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1.  INTRODUCTION 


During  the  analysis  of  120-mm  gun  firings,  designed  to  look  at  the  interior  ballistic  characteristics  of 
combustible  cartridge  cases  (Robbins,  Koszoru,  and  Minor  1986),  XNOVAKTC  (XKTC)  (Gough  1980) 
was  noted  to  be  in  agreement  with  measured  pressure-time  curves  as  well  as  with  pressure-difference  curves, 
while  IBHVG2  (Anderson  and  Fickie,  1987)  gave  a  calculated  maximum  breech  pressure  which  was  42 
MPa  higher  than  measured.  Parametric  studies  were  performed  using  XKTC  to  attempt  to  attribute  this 
disparity  to  the  various  processes  omitted  from  IBHVG2.  The  boattail  instrusion  was  calculated  to  account 
for  14  MPa,  with  effects  of  flamespreading  and  intergranular  stress  accounting  for  3  MPa  each.  Subsequent 
calculations  (Robbins  1986)  indicated  chambrage,  Dropellant  packaging,  wave  dynamics,  and  multiphase 
effects  (the  solid  propellant  velocity  lag  and  concomitant  formation  of  an  ullage  region  between  the  projectile 
base  and  the  propellant  bed)  as  contributors  to  the  differences  between  the  lumped-parameter  and  two-phase 
interior  ballistic  codes.  In  this  report,  we  show  that  the  influence  of  chambrage  and  propellant  velocity  lag 
on  the  pressure  gradient  may  be  represented  in  analytical  form.  We  compare  the  analytical  pressure  gradient 
with  that  predicted  by  XKTC  and  assess  the  extent  to  which  the  effects  of  chambrage  and  propellant  velocity 
lag  account  for  the  differences  in  ballistic  predictions. 

The  NOVA  codes,  of  which  XKTC  is  the  latest  version,  have  been  used  with  uncompromised  data 
bases  to  model  gun  systems  with  much  success  (Robbins,  Koszoni,  and  Minor  1986;  Robbins  1983;  Robbins 
and  Horst  1984).  Since  XKTC  calculates  the  pressure  gradient  from  first  principles  and  agrees  with  gun 
firings,  XKTC  is  assumed  correct.  Accordingly,  all  the  lumped-parameter  computer  runs,  with  different 
gradient  equations,  are  compared  with  equivalent  XKTC  computer  runs. 

1.1  Models.  Gradient  equations  have  been  developed  by  Gough  (Gough,  no  date)  to  look  at  (1)  the 
effects  of  chambrage,  (2)  the  effects  of  propellant  velocity  lag,  and  (3)  both  effects  at  the  same  time.  The 
governing  partial  differential  equations,  assumptions,  and  definitions  are  given  for  each  effect  followed  by 
the  resultant  gradient  equations  and  an  equation  for  the  total  kinetic  energy  of  the  propellant  and  gas.  The 
details  of  the  solution  technique  and  mathematical  procedure  are  documented  in  Appendices  1-3. 

1 .2  Influence  of  Chambrage.  For  the  chambrage  gradient  equation,  the  propellant  is  assumed  to  be 
uniformly  distributed  between  the  breech  and  the  base  of  the  projectile,  and  the  variation  in  area  is  assumed 
to  be  confined  to  the  chamber.  The  basic  results  presented  in  this  section  have  been  previously  derived 
(Vinti  1942)  and,  more  recently,  in  a  somewhat  different  context,  (Morrison  and  Wren  1988). 

The  continuity  and  momentum  equations  for  unsteady  flow  of  a  homogeneous  inviscid  substance  are 


3p  2.  dpAu 

dt  A  dz 
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(C.  01) 
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Making  the  Lagrange  assumption 
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results  in  the  pressure  distribution 
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Substituting  the  projectile  acceleration  given  by 

into  (C.OS)  results  in  a  pressure  distribution  in  terms  of  the  pressure  at  the  base  of  the  projectile  PB 


(C.10) 


p(z)  =  p( 0)  +  a,(r)/,(z)  +  aJitV^a  +  b{t)J2{z), 
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With  the  mean  pressure  Pm  defined  as 


P„  = 


rzr 

P(z)A{z)dz 

rzf 

J  A(z)dz 


(C.14) 


and  using  (C.l  1)  and  identifying  P(0)  with  PBr,  then  the  base  and  breech  pressures  are  given  in  terms  of  the 
mean  pressure  by 


Pm  +  aWfr,)  +  b(ty2(zp) 


aMh (zP) 
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and  . 

PBr  =  pb  ~  Wi(z,>  -  b(t)J2(zp),  (C.16) 


where 

rzr 

h(zP)  =  Jo  Aizy^zyiz 
and 

rz> 

U*p)  =  l  A(z)J2(z)dz. 


The  kinetic  energy  of  the  gas/solid  mixture  is 

K£  m  AlV2pCUzp) 
2 g0V\zp)  ' 


(C.  17) 


(C  .18) 


(C.  19) 


1-3  Influence  of  Propellant  Velocity  Lag.  For  the  two-phase  gradient  equation  developed  in  this 
section,  the  chamber  and  tube  are  assumed  to  have  the  same,  uniform  diameter.  The  coupling  of  the  two-phase 
effects  with  chambrage  is  addressed  in  the  next  section.  The  propellant  is  assumed  to  be  uniformly  distributed 
between  the  breech  face  and  the  leading  edge  of  the  bed  with  the  leading  edge  of  the  bed  being  initially  at 
the  base  of  the  projectile.  However,  once  the  projectile  begins  to  move,  the  leading  edge  of  the  propellant 
bed  is  not  assumed  to  remain  in  contact  with  the  base  of  the  projectile.  Instead,  the  motion  of  the  leading 
edge  of  the  bed  is  explicitly  modeled  by  reference  to  the  equation  of  motion  of  the  solid  phase  which  includes 
the  influence  of  both  pressure  gradient  (buoyancy  force)  and  interphase  drag.  The  existence  of  a  region  of 


3 


ullage  between  the  propellant  bed  and  the  base  of  the  projectile  is  also  recognized  explicitly.  The  requirement 
that  the  gas-phase  density  be  everywhere  uniform  causes  the  gas  velocity  at  the  leading  edge  of  the  propellant 
bed  (Ug)  to  have  a  much  more  complex  time  dependence  than  that  in  the  simple  Lagrange  analysis.  Instead 
of  being  a  fixed  fraction  of  the  projectile  velocity  at  all  times,  the  value  of  Ug  may  be  significantly  greater 
than  the  projectile  velocity  during  the  pressurization  phase  since  the  ullage  requires  a  net  compression  to 
maintain  equilibrium  with  the  conditions  in  the  mixture  region  where  combustion  is  occurring.  The  complex 
time  dependence  of  Ug  will  be  seen  to  have  important  repercussions  in  respect  to  the  behavior  of  the  pressure 
gradient.  It  is  believed  that  this  is  the  first  time  that  these  effects  have  been  modeled  analytically. 


The  continuity  equations  for  the  mixture  region  are 


3ep  depu 
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de  3[(l-e)p  pup] 
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(L.  1) 

(L.2) 


Making  assumptions  analogous  to  the  Lagrange  assumption 
3d 

■J~  =  0  (gas  density  constant  throughout  the  tube)  (L3) 

oz 

3-  =  0  (porosity  constant  throughout  the  mixture  region)  (LA) 


and  noting  the  assumption  m(t)  is  only  a  function  of  time,  then 
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Zb 


and 
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(L.5) 


(L.  6) 


At  the  internal  boundary  defined  by  the  leading  edge  of  the  propellant  bed  we  may  write  the 
macroscopic  balance  of  mass  in  the  following  form 

e  (Us-Up)  =  Us>  -  U„. 

This  result  includes  the  assumption  that  the  density  jump  across  the  boundary  is  negligible,  as  will  always 
be  the  case  if  the  Mach  number  is  small  compared  with  unity,  and  which  is  in  any  case  consistent  with  the 
present  assumption  L.3. 
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Consider  the  momentum  equation  for  the  mixture  region, 
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where  dP^dt  is  determined  numerically,  or  from 
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For  the  ullage  region,  the  momentum  equation  for  the  gas  is 


dP 

dz 


P. 

So 


(du 

u 


+  u 


du} 

dz 


(L.  23) 


and,  therefore,  the  pressure  in  the  ullage  region  is 
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The  kinetic  energy  of  the  propellant  and  gases  is 


KE 


6*7 


W2 


p KL 
bgo 


The  mean  pressure  over  the  mixture  region  is 


1.4  Combined  Influence  of  Cha.r»brage  and  Propellant  VelocitvLag.  The  Robbins-Gough-Anderson 
(RGA)  gradient  equation,  which  combines  the  influences  of  chambrage  and  multiphase  flow,  assumes  that 
the  variation  in  area  is  confined  to  the  propellant  chamber  and  that  the  propeiiant  initially  nils  the  chamber. 
The  cross-sectional  area  within  the  ullage  is  accordingly  uniform  and  equal  to  AB  at  all  points.  Nomenclature 
is  as  in  the  previous  gradient  equations. 

The  continuity  equations  for  the  mixture  region  are 


3ep  3epuA 

dt  +  A  dz 


(/U) 


de  _  1_  3(l-e)ppMpA 

Pp  dt  A  dz 


m(t). 


(R.2) 


Making  assumptions  analogous  to  the  Lagrange  assumption,  namely 


dp  de 

dz  dz 


=  0  in  the  mixture  region 


and 


dz 


=  0  in  the  ullage  region 


(R.  3) 


{R. 4) 


and  noting  the  assumption  m(r)  is  only  a  function  of  time,  then 


and 


u.  = 


u  = 


UpAsVlz) 

V(zb)A(z) 

u^Asm 

V(zb)A{zy 


(R.  5) 


(R.  6) 


At  the  internal  boundary  defined  by  the  leading  edge  of  the  propellant  bed  we  may  write  the 
macroscopic  balance  of  mass  in  the  following  form 

e  (Ue-Up)  =  Ug+-Up 

This  result  includes  the  assumption  that  the  density  jump  across  the  boundary  is  negligible,  as  will  always 
be  the  case  if  the  Mach  number  is  small  compared  with  unity,  and  which  is  in  any  case  consistent  with  the 
present  assumption  R.3. 
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Consider  the  momentum  equations  for  the  mixture  region. 


z[*T  +  +  £S*  I  *  * 


+  rhuB 


i  ram-aw  am-awgl  +  „-c)8. 

A  [  dt  dz  ]  dz 


=  f,~mup, 


where 


and  for 


f,  =  ^77-^  p  f» 

up 


REN  = 


p dM~M 


and  for 


and  with  e*  the  initial  porosity 


0.5  +  Sf 

(i.5  sr. 


Defining 


/«  = 


2.5k  (1-e)  O045  f  p>E 

RENMl  [(1-0  ej 


<J>,  —  <j>  —  p AgLIC 


(R.  10) 


d2lnp  ^  &Abpb  goAfo* 

IF  "  A  )  vfmp  m„vf  ’ 


(/Ml) 
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where 


and 


c,«) 


m 


•  1 

1 

n 

P'V'l 

+ 

v} 

ml 

-v) 

c 

m 

—  + 
Pp 

Pp 

dx 

Pp 

s 

dt’ 

dx 

d2x 

dx 

+ 

P, 

s  y 

(R.  12) 


(*.13) 


or 


with 


and 


dS  (dx)2  c  d2x 

m  =  p'  S  VS)  +  p'  s  a1 


dx  n, 
—  =  <jp; 


A  _  p.-i  dPm 

=  MP”  IT 


(*.  14) 


(K.15) 


where  dPm/dt  is  determined  numerically,  or  from 


ifc/rr  mrT 

dPm  ~  -  KV  +  1=* 


dt 


and 


Then 


f»(z)  =  F*  +  (at(0  +  a&)PB)  Jx{z)  +  &(r)/2(z) 
in  the  mixture  region,  where 


a,(0  = 


CM, 

S.V2(z*) 


<M 

V{z>) 


Vi  +  D  - 


f ffA  „ 

tuf  r  rv 


1- 


p/&*> 


(*.  16) 

(K.17) 


(/?.18) 
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and 


and 


(hit) 


ce^aIPb 

V2(z„)eMp 


1 


fyP_ 

P PV(zb) 


bit ) 


C<Ml 

2g0V\h)' 


J(z)  -  ^ 
2  '  A\zj 


v*  ?  - »; 


dlnp 


dt 


2AB<S>,U.  o 

+  ~VM  +  ♦>  ^  /- 


and 


4>i  -  1  -  4>  -  <(>. 

<t>3  =  4>.  Ul  +  (1-<|»)  U2f, 

Z,#. 

D  =  ~  C,(0, 


£  =  1  - 


LAg 

V,  • 


Also  in  the  ullage  region 


a/» 

dz 


£.  r* 

go  L  9r 


+ 


u 


du' 

dz 


iR.  19) 

(£.20) 

(£.21) 

(£.2  la) 

(£.22) 

(£.23) 


12 


P(Z)  =  P(Z>)  -  Jg  +  *, A)  (z  -2b) 


3/ 


where 


d2lnp  f  dlnf) 


Defining  the  mean  pressure  as 


P_  = 


f  P(z)A(z)dz 
J  A(z)dz 


and  substituting  for  Vp  and  d2  In  p/d/2, 
then 


Par  =  /*.  l^l-^Vl(Zfc)  + 

.  r<««p 


pLAB 

pL'Aj 

Mp 

2VFMp 

l2 

f 

gXPm) 

T* 

m,vf 

-  a,(/y,(zt)  ~  *>{t)J2(zb)  - 


ABpLP n 


(R.  24) 


(7?  .25) 


(R.26) 


p‘  '  ['■  -  £  b>  -  (irj 

_  L„,  ^  sAlP„  _ 

3j.Vfe,)  l  ‘l  '  M,Vr 


ScAbP ru 

MpvF 


ai(t)Ji(zb) 


bm(zb)  AMtVfa)  ABLb(t)J2{zb)  p AlL2Pro 

+  V(zp)  +  V(zp)  +  V(zp)  +  2V{zp)Mp 


-  aWfa)  -  b(t)J2(zb) 


ABpLP n 


ABL02{t)J x{zb)  pAjL2  QlitV&b)  ,  s  r  /  \ 

1  +  — vw —  wm  +  ■w  ~  “,(')/'(2,) 


^  p LAB  p L?Ab  ^  p AbL? a\ 


2  VfMp  3  V(zp)VfMp 


(R.27) 
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where 


Hh)  =  JT  A(zVx(z)dz 


and 


Uh)  =  [ 


'*  V2(z) 


dz. 


The  kinetic  energy  is 

££  = 


A2(z) 

'l-e'j  U2pA2BppUzb)  _e_  tfofp/4(z») 
,  J  V(z>)2  +  2g„  V(z,)2 


p AbL 

"6iT 


3V2  +  3V,p£ 


dlnp 

dt 


+  L 


dlnpy 

dt  ) 


a\ 


The  mean  pressure  over  the  mixture  region  is 


„  „  (<*i(0  +  HtWa)  ,,  x  b(t)JA(zb) 

Pmx  =  PBr  +  - 777. 77 -  A(*»)  + 


V(z*) 


V(z*) 


(£.28) 


(£.29) 


(£.30) 


(£.31) 


2.  CALCULATIONS 


A  number  of  propellant  charges  were  simulated  with  XKTC  to  probe  the  influence  of  chambrage  and 
velocity  lag.  The  same  charges  were  simulated  with  IBRGA  (input  description  and  listing  given  in 
Appendices  4  and  5)  using  data  bases  as  consistent  with  those  of  XKTC  as  the  physical  scope  of  the 
lumped-parameter  model  would  permit,  and  using  each  of  the  various  pressure  gradient  models.  The  resultant 
maximum  breech  pressures,  velocities,  and  histories  of  pressure-time  and  mean-pressure-to-base-pressure 
curves  are  compared. 

The  calculations  perfon  ned  with  XKTC  involved  data  bases  with  evenly  distributed  seven-perforated 
propellant  having  an  initial  porosity  of  0.4,  zero  barrel  resistance,  and  with  all  the  propellant  ignited  at  the 
initial  instant.  All  calculations  were  performed  for  a  flat-based  projectile  and  nominal  heat  loss. 


The  parameters  used  in  the  computer  codes  were: 
Large  chamber  simulations: 

Bore  diameter 
Volume 
Travel 

Propellant  mass 
Projectile  mass 


127  mm 
9,832.2  cm3 
4.572  m 
9.8  kg 

2.45, 9.8,  39.2  kg 
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Small  chamber  simulations: 

Bore  diameter  28.65  mm 

Volume  98.322  cm3 

Travel  1.880  m 

Propellant  mass  .098  kg 

Projectile  mass  .0245,  .098,  .392  kg 

Propellant  characteristics: 

Impetus  1,136  J/g 

Covolume  .976  cm3/g 

Gamma  1.23 

Flame  temperature  3,141  K 

Molecular  weight  23.0  g/gmole 

Density  1.66  g/cm3 

Burning  rate  1 . 105 19p10  mm/s  (p  is  in  MPa) 

Three  maximum  breech  pressures  were  studied,  nominally  517  MPa,  345  MPa,  and  172  MPa.  The 
maximum  pressures  were  achieved  by  allowing  the  web  to  vary,  with  the  grain  length  kept  between  two  to 
three  times  the  outer  diameter.  Calculations  were  performed  for  both  bore  diameter  chambers  and  right 
circular  cylindrical  chambers  with  a  larger  diameter  than  the  bore.  Chambrage  (defined  as  the  distance  over 
which  the  chamber  tapers  down  to  the  bore  area)  was  76.2  mm  for  the  large  volume  chamber,  with  the 
chamber  length  being  changed  from  776.22  mm  (straight  chamber)  to  541.02  mm.  For  the  smaller  volume 
chamber,  we  had  25.4  mm  of  chambrage  with  the  chamber  length  being  changed  from  152.4  mm  (straight 
chamber)  to  101.6  mm. 

3.  RESULTS 

3.1  Chambrage.  The  effect  of  chambrage  was  assessed  at  a  nominal  pressure  of  345  MPa  using 
XKTC.  Calculations  were  performed  for  a  straight  chamber  to  determine  grain  dimensions  corresponding 
to  a  maximum  breech  pressure  of  345  MPa.  These  grain  dimensions  were  then  used  with  a  chamber  with 
the  same  volume  but  hrving  chambrage.  The  presence  of  chambrage  was  found  to  reduce  the  maximum 
breech  pressure  substantially.  Finally,  XKTC  was  run  with  chambrage  and  the  dimensions  of  the  grains 
changed  to  restore  the  maximum  breech  pressure  to  345  MPa.  These  sets  of  grain  dimensions  were  then 
used  with  a  version  of  IBRGA  into  which  the  chambrage  gradient  equation  had  been  encoded.  Equivalent 
computer  runs  on  IBRGA  with  the  Lagrange  gradient  are  also  given  in  Table  1. 

Adding  chambrage  to  the  gradient  equation  in  the  lumped-parameter  code  captures  the  large  drop  in 
maximum  breech  pressure  and  the  concomitant  drop  in  muzzle  velocity  seen  when  chambrage  is  added  to 
XKTC.  The  gradient  equation  with  chambrage  also  gives  the  proper  amount  of  increase  in  maximum  breech 
pressure  and  muzzle  velocity,  as  seen  in  XKTC  calculations,  when  the  grain  dimensions  are  changed  to 
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TABLE  1 


Chambrage  Calculations 


Chamber 

volume, 

cm1 

c/m 

Chambrage, 

mm 

XKTC 

Chambrage 

gradient 

Lagrange  gradient 

Maximum 

breech 

pressure, 

MPa 

Muzzle 

velocity, 

m/s 

Maximum 

breech 

pressure, 

MPa 

Muzzle 

velocity, 

m/s 

Maximum 

breech 

pressure, 

MPa 

Muzzle 

velocity, 

m/s 

9,832.2 

0.25 

none 

345 

769 

34  5 

775 

350 

775 

9,832.2 

0.25 

76.2 

336 

764 

337 

766 

9,832.2 

0.25 

76.2 

345 

llH 

346 

776 

358 

784 

9,832.2 

1.0 

none 

345 

1,352 

347 

1,363 

347 

1,363 

9,832.2 

1.0 

76.2 

310 

1,303 

307 

1,308 

9,832.2 

1.0 

76.2 

345 

1,376 

339 

1,375 

384 

1,429 

9,832.2 

4.0 

none 

345 

1,964 

360 

1,962 

360 

1,962 

9,832.2 

4.0 

76.2 

268 

1,764 

273 

1,771 

9,832.2 

4.0 

76.2 

345 

1,964 

339 

1,990 

452 

2,184 

98.322 

0.25 

none 

345 

895 

348 

880 

348 

880 

98.322 

0.25 

25.4 

337 

895 

334 

872 

98.322 

0.25 

25.4 

345 

901 

341 

879 

355 

886 

98.322 

1.0 

none 

345 

1,560 

347 

1,517 

347 

1,517 

98.322 

1.0 

25.4 

308 

1,526 

305 

1,469 

98.322 

1.0 

25.4 

345 

1,595 

340 

1,536 

387 

1,582 

98.322 

4.0 

none 

345 

2,298 

371 

2,198 

371 

2,198 

98.322 

4.0 

25.4 

275 

2,149 

280 

2,028 

98.322 

4.0 

25.4 

345 

2,354 

342 

2,230 

457 

2,389 
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restore  the  maximum  breech  pressure  to  345  MPa.  If  the  Lagrange  gradient  equation  were  to  be  used,  a 
difference  from  XKTC  of  approximately  10%  at  a  charge  mass  to  projectile  mass  ratio  (c/in)  of  one  and 
approximately  30%  at  a  c/m  of  four  would  occur  in  the  predicted  values  of  maximum  breech  pressure. 


The  chambrage  gradient  equation  gives  calculated  maximum  breech  pressures,  which  are  as  much  as 
7%  high  and  muzzle  velocities  low  by  3  to  4%  when  compared  to  XKTC  calculations.  The  resultant  gradient 
structure  is  illustrated,  in  Figures  1-6,  with  representative  plots  of  the  ratio  of  the  mean  pressure  to  the  base 
pressure  for  both  XKTC  and  the  analytic  chambrage  gradient  equation.  The  figures  also  include  histories 
of  breech  and  base  pressure. 

3.2  Multiphase  Flow.  The  effects  of  multiphase  flow  were  examined  using  XKTC.  Calculations 
were  based  on  a  straight  chamber,  nominal  maximum  breech  pressures  of  172  MPa,  345  MPa,  and  5 17  MPa; 
large  and  small  chamber  volumes;  and  values  of  c/m  equal  to  .25, 1,  and  4,  The  grain  dimensions  from  the 
XKTC  calculations  were  used  to  define  equivalent  data  bases  for  IBRGA  into  which  the  two-phase  gradient 
equation  had  been  encoded. 

The  pressure  at  which  the  propellant  bums  is  taken  as  the  mean  pressure  over  the  mixture  region 
instead  of  the  mean  pressure  over  the  entire  volume  behind  the  projectile,  as  is  conventionally  assumed. 
Comparable  Lagrange  gradient  calculations  are  also  supplied  in  Table  2. 

The  calculated  maximum  breech  pressures  and  velocities  forboththe  two-phase  and  Lagrange  gradient 
for  c/'m  of  .25  are  in  very  close  agreement  with  XKTC  at  all  pressures.  For  calculations  at  c/m’s  of  1  and 
4,  the  Lagrange  and  two  phase  gradients  tend  to  give  close  maximum  breech  pressures,  at  worst  about  6% 
different  from  XKTC. 

The  history  of  the  ratio  of  space  mean  to  base  pressures  for  representative  two-phase  gradient 
calculations  and  for  corresponding  XKTC  calculations  is  plotted  in  Figures  7-12,  together  with  associated 
breech-pressure  and  base-pressure  histories.  The  Lagrange  gradient  would  be  a  straight  line  with  a  ratio 
value  of  1.083  for  a  c/m  of  0.25;  1.333  for  a  c/m  of  1;  and  2.333  for  a  c/m  of  4.  The  two-phase  gradient 
equation  has  the  same  shape  and  magnitude  of  the  first  peak  in  the  mean  pressure  to  base  pressure  ratio  as 
XKTC,  with  the  maximum  occurring  earlier  in  the  analytic  two-phase  gradient  calculations  and  the  recovery 
to  the  Lagrange  ratio  value  not  being  as  fast  as  XKTC.  The  comparison  of  shape  and  magnitude  of  the 
pressure  ratio  histories  strongly  suggests  that  the  dominant  physical  processes  are  captured  by  the  two-phase 
gradient  equation.  In  particular,  the  undulatory  character  of  the  curve  is  seen  to  be  attributable,  at  least  in 
part,  to  the  two-phase  aspects  of  the  flow,  rather  than  being  due  to  transient  or  wave  propagation  phenomena 
as  might  at  first  be  thought. 
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P-Ratio  Mean/Base  ^  P-Ratio  Mean/Base 


Pressure  (MPa)  Pressure  (MPa)  Pressure  (MPa) 


TABLE  2.  Two-Phase  Calculations 


Chamber 

volume, 

cm1 

c/m 

XKTC 

Two-Phase  gradient 

Lagrange  gradient 

Maximum 

breech 

pressure, 

MPa 

Muzzle 

velocity, 

m/s 

Maximum 

breech 

pressure, 

MPa 

Muzzle 

velocity, 

m/s 

Maximum 

breech 

pressure 

MPa 

Muzzle 
,  velocity, 
m/s 

9,832.2 

0.25 

174 

549 

177 

557 

179 

557 

9,832.2 

0.25 

345 

769 

348 

773 

350 

775 

9,832.2 

0.25 

518 

887 

521 

891 

519 

891 

9,832.2 

1.0 

173 

957 

177 

970 

184 

966 

9,832.2 

1.0 

345 

1,352 

343 

1,355 

347 

1,363 

9,832.2 

1.0 

516 

1,593 

498 

1,556 

497 

1,560 

9,832.2 

4.0 

172 

1,389 

172 

1,380 

184 

1,354 

9,832.2 

4.0 

345 

1,964 

351 

1,930 

360 

1,962 

9,823.2 

4.0 

518 

2,391 

524 

2,290 

528 

2,306 

98.322 

0.25 

172 

641 

177 

632 

179 

630 

98.322 

0.25 

345 

895 

346 

880 

348 

880 

98.322 

0.25 

518 

987 

514 

974 

513 

974 

98.322 

1.0 

172 

1,130 

180 

1,107 

186 

1,094 

98.322 

1.0 

345 

1,560 

343 

1,516 

347 

1,517 

98.322 

1,0 

516 

1,756 

510 

1,703 

512 

1,705 

98.322 

4.0 

172 

1,682 

174 

1,581 

185 

1,524 

98.322 

4.0 

345 

2,298 

360 

2,188 

371 

2,198 

98.322 

4.0 

518 

2,685 

545 

2,516 

552 

2,527 
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Lgure  7 y  XKTC  Calculation  With 
Large  Volume ,  c/m  of 
0.25,  345  MPa  and  No 
Chambraqe . 


igure  8 XKTC  Calculation  With 
Large  Volume,  c/m  of 
1.0,  345  MPa  and  No 
Chambraqe. 
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Figure  9.  XKTC  Calculation  With 
Large  Volume,  c/m  of 
4.0,  345  MPa  and  No 
Chambrage. 


Figure  10.  IBRGA  Calculation  With 
Two-Phase  Gradient 
With  Figure  7  Values. 
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Figure  11.  IBRGA  Calculation  With 
Two-Phase  Gradient 
With  Figure  8  Values  * 
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Figure  12.  IBRGA  Calculation  With 
Two-Phase  Gradient 
With  Figure  9  Values., 
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The  drop  off  of  the  pressure  ratio  at  the  slivering  event  seen  in  IBRG  A  calculations  may  cause  numerical 
problems.  (This  event  is  smoothed  over  in  XKTC  since  all  the  grains  do  not  sliver  at  the  same  time  because 
of  the  different  pressure  fields  seen  by  the  grains,  as  well  as  the  finite  time  required  for  the  transfer  of 
information.)  This  problem  is  not  felt  to  affect  the  calculated  velocity  significantly,  but  a  smoothing  routine 
has  been  incorporated  into  IBRGA  to  circumvent  any  difficulty.. 

The  two-phase  gradient  equation  also  has,  as  input,  a  multiplier  to  the  friction  ft  ctor.  By  making  this 
value  small,  an  approximation  to  stick  propellant  behavior  can  be  captured. 

3.3  Combined.  Table  1  is  reproduced  in  Table  3  with  the  chambrage  gradient  equation  replaced  by 
the  RG  A  gradient  equation  with  the  appropriate  kinetic  energy  term.  As  with  the  two-phase  gradient  equation, 
the  burning  rate  is  calculated  from  mean  pressure  over  the  mixture  region. 

The  shapes  of  the  ratio  and  pressure-time  plots,  Figures  13-18,  are  qualitatively  similar  to  those 
produced  by  XKTC.  This  reinforces  the  idea  that  both  two  phase  and  chambrage  effects  are  important  factors 
in  respect  to  the  time-dependent  behavior  of  the  pressure  gradient.  The  result  of  the  RGA  gradient  equation 
calculations  are  in  general  slightly  better,  when  compared  to  XKTC,  than  the  corresponding  calculations 
with  the  chambrage  gradient  equation. 

One  of  the  RGA  gradient  equation  inputs  is  a  multiplier  to  the  friction  factor.  By  making  this  value 
small,  an  approximation  to  stick  propellant  behavior  can  be  captured. 

4.  CONCLUSIONS 

It  is  concluded  that  incorporating  chambrage  into  a  gradient  equation  is  important  and  can  make  as 
much  as  a  20  to  30%  difference  in  predicted  maximum  breech  pressure  at  c/m  of  4.  The  inclusion  of  two-phase 
effects  into  the  gradient  equation  seems  to  be  well  motivated  in  that  a  large  portion  of  the  structure  of  the 
gradient  history  is  captured. 

At  low  c/m’s  (0.25),  any  of  the  gradient  equations  will  give  good  results;  but  as  the  c/m  gets  larger 
(>1),  the  correction  for  chambrage  is  required. 

The  use  of  the  RGA  gradient  equation  offers  a  natural  way  to  account  for  stick  propellant  phenomena 
having  to  do  with  stick  propellant  motion  and  the  pressure  history  in  which  it  bums. 

The  use  of  the  proper  kinetic  energy  term  for  the  gradient  equation,  as  well  as  the  calculation  of  the 
mean  pressure  over  the  mixture  region,  should  be  included  in  any  gradient  model  when  applicable. 
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TABLE  3.  RGA  Calculations. 


Chamber 

volume, 

cm3 

Chambrage, 
c/m  mm 

XKTC 

RGA  gradient 

Lagrange  gradient 

Maximum 

breech 

pressure, 

MPa 

Muzzle 

velocity, 

m/s 

Maximum 

breech 

pressure, 

MPa 

Muzzle 

velocity, 

m/s 

Maximum 

breech 

pressure, 

MPa 

Muzzle 

velocity, 

m/s 

9,832.2 

0.25 

none 

34  5 

769 

348 

773 

350 

775 

9,832.2 

0.25 

76.2 

336 

764 

336 

766 

9,832.2 

0.25 

76.2 

345 

774 

345 

775 

358 

784 

9,832.2 

1.0 

none 

345 

1,352 

343 

1,355 

347 

1,363 

9,832.2 

1,0 

76.2 

310 

1,303 

304 

1,300 

9,832.2 

1.0 

76.2 

345 

1,376 

337 

1,369 

384 

1,429 

9,832.2 

4.0 

none 

345 

1,964 

351 

1,930 

360 

j.961 

9,832.2 

4.0 

76.2 

268 

1,764 

266 

1,747 

9,832.2 

4.0 

76.2 

345 

1,964 

334 

1,964 

4  52 

2,183 

98.322 

0.25 

none 

345 

895 

346 

879 

348 

880 

98.322 

0.25 

25.4 

337 

895 

332 

872 

98.322 

0.25 

25,4 

345 

901 

340 

878 

355 

886 

98.322 

1.0 

none 

345 

1,560 

342 

1,515 

347 

1,517 

98.322 

1.0 

25.4 

308 

1,526 

302 

1,468 

98.322 

1.0 

25.4 

345 

1,595 

337 

1,535 

387 

1,582 

98.322 

4.0 

none 

345 

2,298 

361 

2, 190 

371 

2,198 

98.322 

4  -  0 

25.4 

275 

2,149 

272 

2,023 

98.322 

4.0 

25.4 

345 

2,354 

336 

2,222 

457 

2,389 
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P-Ratio  Mean/Base  _  P-Ratio  Mean/Base 


Pressure  (MPa)  Pressure  (MPa)  Pressure  (MPa) 


Intentionally  lot  blank. 
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Intentionally  left  blank. 
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6.  LIST  OF  SYMBOLS 


e,  porosity  defined  as  the  volume  of  gasAotal  volume 
X,  particle  shape  factor  for  gas  flow  resistance 
p,  viscosity  of  gas 

p,  density 

pp,  propellant  density 

mass  fraction  burned 
A(z),  area  at  position  z 

Ab,  tube  area 

C,  total  charge  mass 

Dp,  effective  diameter  of  the  propellant  grain,  i.e„  6*  volume  of  propellant  grain/surface  area  of 
a  propellant  grain 

f$,  interphase  drag 

/so,  baseline  friction  factor 

/se,  dimensionless  friction  factor,  a  function  of  Reynolds  number  and  particle  shape 

GD,  grain  diameter 

GL,  grain  length 

g0,  a  constant  to  reconcile  dimensions 

L,  the  distance  between  the  leading  edge  of  the  propellant  bed  and  the  base  of  the  projectile, 
i.e.,  L  =  (Zp-Zfc) 

m(t),  rate  of  combustion  of  propellant 
m,  mass  of  propellant  burned 

M, ,.  mass  of  projectile 

mw,  molecular  weight  of  propellant  gas 

P,  pressure 

PB,  projectile  base  pressure 

PBr,  breech  pressure 

Pm,  mean  pressure 

Pra,  resistive  pressure 

R,  universal  gas  constant 
REN,  particle  reynolds  number 

S,  surface  area  of  the  propellant 

T,  gas  temperature 

t,  time 

u,  velocity 
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U|t  gas  velocity  on  the  mixture  side  at  the  leading  edge  of  the  propellant  bed,  i.e.,  u^) 

U,„  gas  velocity  on  the  ullage  side  of  the  leading  edge  of  the  propellant  bed 

up,  propellant  velocity 

Up,  velocity  of  the  leading  edge  of  the  propellant  bed,  i.e.,  up(Zt,) 

V(z),  volume  between  the  breech  face  and  position  z 

VF,  free  volume  for  gases  to  expand  into 

VOI  initial  empty  volume  of  chamber 

Vp,  velocity  of  the  projectile 

x,  distance  burned  into  the  propellant 

z,  axial  position  measured  from  the  breech  face 

Zp,  position  of  the  leading  edge  of  the  propellant  bed 

Zp,  position  of  the  base  of  the  projectile 

zp#,  initial  position  of  the  base  of  the  projectile 

•,  derivative  with  respect  to  time 

second  derivative  with  respect  to  time 
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APPENDIX  1: 

INFLULNCE  OF  CHAMBRAGE 
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* 


Intentionally  left  blank. 
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The  following,  is  work  by  Gough  (Gough,  in  preparation)  with  modifications  and  elaborations  by 
Robbins. 

Influence  of  Chambrage.  For  the  effects  of  chambrage  on  the  form  of  the  gradient  equation,  the 
propellant  is  assumed  to  be  uniformly  distributed  between  the  breech  and  the  base  of  the  projectile,  and  all 
variations  in  tube  area  arc  confined  to  the  chamber. 

The  continuity  and  momentum  equations  for  unsteady  flow  of  a  homogeneous  inviscid  substance 
through  a  tube  with  variable  area  are 


§P 

dr 

du 

dt 


1  fdpAu^ 

A  {  dz  J  ~ 


du 

+  <>“&  +  *• 


a p 

dz 


=  0. 


Making  the  Lagrange  assumption 


3p 

dz 


=  0, 


then  C.l  becomes 


dp 

dr 


P 

A 


(t) 


The  density  may  be  written  as 


P  = 


V(zfY 


Differentiating  (C.5),  we  get 


dp 

dt 


C  )  (dV£)) 

l^,)J  V  dt  j 


But, 

and 


V(zp)  = 


dV(zp) 

dt 


V(2po)  + 

=  ABip  =  AgVp 


Ab 


(C.l) 

(C.2) 

(C.3) 

(C.4) 

(C.5) 


(C.6) 


(C.7) 

(C.8) 


31 


therefore, 


But, 


dp  CAbVp 

to  ”  V\z„ ) ' 


P 


C 

v(zPy 


(C.9) 


therefore, 


dp  ^»P  Vp 

to  "  V(z,)  • 


Substituting  (C.10)  into  (C.4),  we  get 


ABpVP  _  p  ( <Mu 
V(zr)  A  V  dz  J 


ABAVp  ~  d Au 
V(zp)  ~  dz  - 


(C.10) 


(C  11) 


(C.12) 


Integrating  (C.12)  with  u(0)  =  0  and  noting  that  Vp  and  V  (Zp)  are  functions  of  time  only. 


A(z)u(z)  =  Au  =  f  A(z)dz  = 

K Jo 


AbVpV(z) 
V(zp)  ■ 


(C.  13) 


Therefore,  the  velocity  u  at  a  given  time  is  given  by: 


Note  that  if  A  =  AB,  then  V(z)  =  A„z,  then  u  =  (z/Zp)  Vp  as  for  the  Lagrange  gradient 
The  pressure  distribution  follows  from  the  momentum  equation  (C.2) 

dP  fdu  du ) 

80  dz  ~  ~  P  [dt  +  Udz} 


(C.14) 

(C  .15) 


(C.16) 


Differentiating  (C.14),  we  get  duldt  and  duldz  for  substitution  into  (C.16). 


To  get  duldt  we  diffetentiate  (C.14),  noting  that  only  Vp  and  V(Zp)  are  functions  of  time  and 


dV{z„)  3(V(z.J  +  (*,-*„>  AB)  ,  „ 

dt  ~  dt  ~  *  p  “  B  p’ 


(CM) 


then 


du 

dt 


ABVpV(z) 

A(z)V(zp) 


abvpv(z)  ^  rav^ 

^a(z)v2{zp) J  I,  a t  ) 


AbV„V(z)  _  A2BVjV(z) 
A(z)V(zp)  A(z)V\zp) 


(C.  18) 


To  get  a uldz,  we  differentiate  (C.14),  noting  that  only  V(z)  and  A(z)  arc  functions  of  z  and 


av(z) 

dz 


)j  A(z)dz 


dz 


=  A(z), 


then 


du 

dz 


AbVpM z) 
V(zp)A(z) 

A,\ r, 
V(zp) 


A„VpV(z) 
A2(z)V(zp)  J 

f  ABVpV(z) > 
A\z)V(zp). 


(« 


m 


Substituting  (C.5),  (C.14),  (C.18),  and  (C.20)  into  (C.16)  yields 


a p 

dz 


CA„V(z)Vp  CV\z)VlAl^ 
A(z)V\zp)  +  A\z)V\zp) 


Integrating  (C.21),  and  noting  Vp  and  V(Zp)  arc  functions  only  of  time,  then 

go  [  &  =  Sow  -  pm 

cabK  r  v(z^ 


Vp  r 

V\zp)  Jo  A(z) 

CV^Al  r*  (V^)\  ('aA') 
\zp)  Jo  U3(Z)J  U2J 


dz. 


(C.  19) 


(C.20) 


(C.21) 


(C.22) 
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Noting 


r  ( 


V\Z)\  f*A\ 


A\z) 


dA 

\dz  j 


dz 


r  m 

i  A(z) 


dz  -  y2(2) 

A{z)  2 A\z) 


(from  integration  by  parts),  then  (C.22),  the  pressure  distribution  behind  the  projectile  is 


8.  (P(z)  -  P(0))  = 


CAjV2  CAbVp 

{ 


V\zp) 


p>  j 


cv\a2b 

2V\zp) 


AOO> 


where 


m  =  £ 


*  m 

A(z) 


dz, 


which  can  be  evaluated  algebraically  or  numerically  if  V(z)  and  A(z)  are  known,  and 


r  y2(z) 


Rewriting  (C.23)  in  the  form 


P(z)  =  P{ 0)  +  a  (07,(2)  +  b(tV2{z) 


gives 


and 


°w  ■  i" 


(CA2BV2p  _  CAs 9, 
v  V%)  V\zp) 


•ZAZ 


6(0  =  - 


cvfa 

2g,V3(zp/ 


Defining  the  mean  pressure  in  the  chamber  to  be 


rv 

J  P(z)A(z)iz 
J  A(z)dz 


(C.23) 

(C.24) 

(C.25) 

(C.26) 

(C.27) 

(C.28) 

(C.29) 
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then  using  (C.26), 


or 


where 


and 


P-  = 


£  [P( 0)  +  am<*)  +  bm(z)}  A(z)dz 


P_  = 


I" 

A(z)dz 

P(°) 

p  A(z)dz 

V(zp) 

<■«)  J 

P  J\(z)A(z)dz 

c 

V(zp) 

»(«  J 

H  Uz)A(z)dz 

o 

V(zp) 


bit) 


F-  =  n0)  +  $7)  +  «’-> 


h(z,)  =  £'  A(z)Jx(z)dz 


■  f  (r 


Z(*) 

A(*) 


dx 


A(z)dz 


J*(zp)  =  J['  A(z)J2(z)dz 
*'•  V2(z) 


■r 


A(z) 


dz, 


(C.30) 


(C.31) 


(C.32) 


(C.33) 


(C  .34) 


and  where,  if  A(z)  and  V(z)  are  known,  J3(Zp)  and  J4(Zp)  may  be  evaluated  algebraically  or  numerically.  The 
acceleration  of  the  projectile  is  given  by 


K 


SqAs 


l Pb-U 


(C.35) 


Substitution  into  (C.27)  gives 

a(0  =  0,(1)  + 


(C.36) 
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where 


,,  CA„ 

~  _  t/2. 


goV\zp)  ^  V(zp) 


AbV2p  SoABPr, 
+  M, 


and 


(h(t)  =  - 


CAj 

M,V\zp) 


Substituting  (C.36)  into  (C.26)  and  (C.32)  and  noting  P*  =  P(0), 


then 


and 


P*  =  [\-(h(mzp)]PB-amzp)  -  b(t)J2(zp) 


P  -  P  +  j  (2  )  +  ML  J(Z) 

*Br  +  \/(„  \  \Zp)  +  t//„  \  *M \?p)* 


V(zp) 


V(2p) 


Substituting  (C.39)  into  (C.40)  and  rearranging  terms. 


P_  = 


1  —  ^f)/x{zp)  "*■ 

-  b(t)J2{zp)  + 


a&W?p) 


Pt  -  ax(t)Jx(zp) 


V(zp) 

ax(tMzP)  b(tyjzp) 


V(zp) 


V(zp) 


(C.  37) 


(C.  38) 


(C.39) 

(C.40) 


(C.41) 


Solving  (C.41),  we  getPB  in  terms  of  Pm 


P. 


Pm  +  ax{t)Jx{zp)  +  b(ty2(zp) 


oi (jMp)  _  mr 4(2,) 

V(ZP)  V(zp) 


1  -  <hM lib)  + 


a&V&j) 

V(zf) 


(C.42) 


Therefore,  assuming  the  mean  pressure,  Pa,  can  be  calculated  from  an  equation  of  state,  then  (C.42)  will 
give  the  projectile  base  pressure,  PB.  The  breech  pressure,  PBr,  is  calculated  from  (C.39)  knowing  PB.  The 
pressure  distribution  is  given  by  (C.26). 


The  evaluation  of  Ji(Zp),  ^(ZpX  and  J^Zp)  can  be  simplified  by  noting  that  the  variation  in  area 
is  confined  to  the  chamber  and,  therefore, 


and 


H2,) 


V(z)  =  Viz^)  +  Aa(z-z ),  for  z>z„, 


-  f  Si  *  -  f  33  *  +  fa* 

Jo  A(z)  Jo  Afz)  J,  A(z) 


A(z) 

v(zpo)  +  AB{z-zpo) 


A(^po)  ^  r  .  iu 

•'*„  ab 

=  w  +  it  (v(g  (^-^)  + 


(C.43) 


(C.44) 


-  V(z)L. - n - • 


(C  .45) 


*-■  -  r  (r  g  *)  *» 

■  r  (r  s  *) «.  •  r  (r  %  • 


A„dz 


=  U*„)  +  J"  ^/,(zg  +  J-  ^V(zp 


>)(z  Zpo)  "*■  (Z  ~Zpo) 


W 


)) 


AB  dz 


=  U2„)  +  VlM  (Zp-Zpo)  +  -  g>-  ^  2p°)2 
+  7  (2p-2pcf> 


(C. 46) 


=  f'  I&>  *  _  fV  V?fe)  .  f*  Rzg  +  Aa(z-zM)f 
Jo  A(z)  02  Jo  A(z)  *  +  J, - £ - —  & 


=  UZpo)  +  tV^>  1  Ab(2,-2po)?-V\z„) 


3 Al 


(C  .47) 


Equations  (C.44)  -  (C.47)  require  the  evaluation  of  the  integrals  J^zg  -  J4(Zpo)  only  once. 
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The  kinetic  energy  (KE)  of  the  gas/solid  mixture  will  be  required,  and  since  p  dV  =  dm  and 
dV  =  A(z)  dz,  then 


SoKE 


I  J'  u2pdV  =  \  ['  u2pA(z)dz 


and  from  (C.5), 


P 


C 

v(z,y 


and  from  (C.14) 


AbVpV(z) 

V(zp]A(z)’ 


therefore, 


or 


1  f*'  CAfryjz)  & 

2  Jc  A(z)V\zp) 

cJp-  if  os  4 

2V\zp)  U>  A(z)  J 


KE 


2g0V3(z,) 


(C  .48) 


(C.49) 
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APPENDIX  2: 

INFLUENCE  OF  PROPELLANT  VELOCITY  LAG 
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Intentionally  left  blank. 
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The  following,  is  work  by  Gough  (Gough,  in  preparation)  with  modifications  and  elaborations  by 
Robbins. 

Influence  of  Propellant  Velocity  Lag.  For  the  effects  of  the  propellant  velocity  lag  on  the  form  of 
the  gradient  equation,  the  propellant  is  assumed  to  be  uniformly  distributed  between  the  breech  face  and 
the  leading  edge  of  the  bed,  with  the  leading  edge  of  the  bed  being  initially  at  the  base  of  the  projectile. 
The  tube  diameter  is  the  same  as  the  bore  diameter. 


The  continuity  equations  for  each  of  the  phases  in  the  mixture  region  are 


3ep  3epu  _ 
a<  +  dz  - 


de 
p'  dt 


3  [(1-e)  p pu„] 
dz 


=  rh(t). 


(LA) 
(L.  2) 


Making  assumptions  analogous  to  the  Lagrange  assumption, 

3p 

=  0  (gas  density  constant  throughout  the  tube)  (L. 3) 

3e 

5-  =  0  (porosity  constant  throughout  the  mixture  region )  (LA) 


and  noting  the  assumption  that  m(t)  is  only  a  function  of  time,  then  from  (L.2), 


3e 

dt 


3  (1-e)  pp 
3 z 


.  v  dur 

(1-£)  P"  Hz  =  rh{t)' 


Noting  (L.4),  then 


Therefore, 


3  (1-e)  P, 

3z 


du. 

(1-e) 


3e 

dt 


0. 

m(t)  =  m 


(L.5) 


(L.6) 


(LA) 
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and 


»,(z)  =  = 


P,  g  ~  *(0 
0-e)  P, 


*  +  /,(0. 


M 


Since  Up(0)  =  0,  then  f\(t)  =  0,  and 


PP  *  -  *(0 

uP(2„)  Ur  (j_e)  Pf  *»’ 


or 


Up  m  Pp  »  -  g(0 
(1-e)  p. 


(£-.9) 


Substituting  (L.9)  into  (L.8)  gives 


U. 

u  (z)  =  —  z. 
zb 


From(L.l), 


dcp  f 

dr  +  U  [ 


de 

dz 


P  ST  +  e  +  ep 


3p 

dz 


du 

dz 


=  m(t), 


(£..10) 


and  using  (L.3)  and  (L.4),  then 


d u 
dz 


T  +  *(0 
ep 


(< a  function  of  time  only). 


Integrating  (L.11)  gives 


=  Jj 


u(z)  =  1 9a  = 


ep 


2  +  /2(0- 


(£..11) 


(£,.12) 


Since  u(0)  =  0,  then  f2(t)  =  0,  and 


or 


u(zb)  -  U,  = 


Hi 

Zb 


~  fr  +  rMf) 
ep 


IT  +  M*) 


ep 


(£..13) 
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Substituting  (L.13)  into  (L.12)  gives 

V . 

u(z)  =  —  z._  (L.14) 

At  the  internal  boundary  defined  by  the  leading  edge  of  the  propellant  bed  we  may  write  the 
macroscopic  balance  of  mass  in  the  following  form 


e  (Ut-U„)  =  Ut+-Up. 


(L.  15) 


This  result  includes  the  assumption  that  the  density  jump  across  the  boundary  is  negligible,  as  will  always 
be  the  case  if  the  Mach  number  is  small  compared  with  unity,  and  which  is  in  any  case  consistent  with  the 
present  assumption  L.3. 


Also,  in  the  ullage  region,  the  continuity  equation  is  given  by 


3p  5p u 
dt  +  dz 


Using  (L.3),  we  get 


du 

dz 


l 

P 


f£\ 

Kdt  j 


(a  function  of  time  only). 


Integrating  (L.17)  we  get 


Using  the  boundary  condition 


gives 


u(z) 


1  ( dp'j 

p  {dt  j 


z  +  m. 


u(z„)  =  V, 


“(O  =  vr  +  « 


Evaluating  (L.20)  at  zb  gives 


U  =  V  +  —  [  -  v  +  L 

'  +  p  [dt  J  V>  +  L  dt  ■ 


Eliminating  Ug+  from  (L.15)  with  (L.21)  gives 


dlnp 


e  V'-V,)  =  Vp  +  L  a  ~  u, 


(L.  16) 

(L.17) 

(L.18) 

(L.  19) 

(L.20) 

(L.21) 

(L.  22) 
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or 


U.  = 


V.  +  i  k  +  L  -  »,) 

■fr 


Differentiating  (L.23)  yields 


Og  =  -  (Jp 


1-e 


V,  L 

+  —  +  — 
e  e 


v  e 

f  J2 


_  1  f  y  +  , 

e2  l  p  dt 


-  £/. 


<rlnp 
"  dr2 


+  — 
e 


dlnp 

dt 


(L.23) 


(L.  24) 


We  anticipate  the  need  to  factor  out  the  base  pressure  PB  dependence  of  both  Vp  and  d2lnp/dt2 . 


The  projectile  acceleration  is  given  by 


%<AbPb 

Mp 


(L.  25) 


To  get  the  base  pressure  PB  dependence  associated  with  d2lnp/-ir2,  we  start  with  the  definition  of 
the  gas  density 

p  -  *  «..2« 


where  m  =  mass  of  gas  burned  and  VF  =  free  volume.  Therefore, 


VF  -  Va  +  Ab  (2p-Zpo)  -  ~  +  J  (£-.27) 

Pp  Pp 

and 

K  -  V,  +  £  =  +  £  (£-.28) 

Pp 

and 

V,  =  -4*P,  +  £  (£-.29) 

P/> 

m  ^ 

where  it  is  understood  Jiat  m!pp,mlpp,  and  m/p.  may  be  in  the  form  Xr1,  Z— ,  ,  where  the  i 

•>,  k  ■V. 

index  refers  to  different  propellants. 
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From  (L.26)  we  get 


lnp  =  lnm  - 


(£-.30) 


d/np  _  m 
dt  m 


(£-.31) 


d2lnp 

dt2 


i  + 


(£-.32) 


Substituting  (L.29)  and  (L.25),  we  get 


T-  =  c&) 


where 


d2lnp 

dr2 


c‘(0  =  *  m 


go^B  g-A^re 

VfA/p  +  ’ 


^  m2 


(£-.33) 


We  still  require  m,  which  we  get  from 


m  =  pp  S 


(£-.34) 


where  S  is  the  surface  area  of  the  propellant  and  dx/dt  is  the  burning  rate  of  the  propellant,  where  dx/dt  is 
of  the  form 


dx 

—  =  aP" 
dt 


Therefore,  from  (L.34), 


A  dx  _  d2x 

m  -  *  s  *  +  *  s  s? 


(L.35) 


dS  (dx Y  „  d2x 

P'  A  U  +  p-  s  aF’ 


and  from  (L.35), 


dr2  dr 


(L.36) 


(£-.37) 
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where  dPJdt  is  determined  numerically  or  from 


ihKT 

dPn  m, 

dt 


-  PmV  + 


mPt 


Also  defining  <|>  as  the  mass  fraction  of  propellant  burned,  then 


P  = 


($>C  —  pABL  <$>.C 


eV(zfc)  eV(z») 


where 


and 


<!>•  =  $  - 


p AgL 


p V  = 


d-»)C 


(1  -e)V(z») 


Differentiating  (L.39),  we  get 


4.  = 


♦  -  £  I L  *  L 


d/rtp'j 

*  / 


Solving  (LAO)  for  e,  we  get 


e  =  1  - 


OHOC 

PfV(zh) 


Differentiating  (L.42)  and  noting  that 


^(Zk)  ~  ^(Zpo)  +  (Zfc  ~  Zpo)  Ag 


and,  therefore. 


then 


^(Zb)  ~  ^&Ag  ~ 


UPA„ 


e  = 


+  (i-«  CV(zb) 


P,V{zb)  p  pV\zb) 


«j)C  (1-fl  Ct/y4a 
PpP^)  +  ppV2(z4) 


(£-.38) 

-39) 

(£.40) 

(£.41) 

(£.42) 

(£.43) 

(£.44) 

(£.45) 
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We  now  have  the  quantities  which  are  required  for  the  solution  of  the  momentum  equation  from  which 
we  get  the  pressure  distribution  in  the  mixture  region. 

Consider  the  momentum  equations  for  each  phase 


where 


and 


where 


and 


(£.46) 

(LAI) 

(£48) 

(£.49) 

(£.50) 


Adding  (L.46)  and  (L.47),  we  get 


or 


8epu  8  (1-e)  p pup  8  epir2  8  (1-e)  p puj  /r 

8r  8r  8z  '  8z  (  ) 

^  [epu  +  (1-e)  p pup]  -  £  [ep"2  +  (1_e)  (£-52) 
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From  (L.10)  we  get  up(z)  =  Up[z/zb],  from  (L.14)  we  get  u(z)  =  U,[z/zb],  from  (L.38)  we  get 


ep 


p 'AgLt 


and  from  (L.40)  we  get 


O-eJP,  = 


V(zb)  V(zb)’ 

(1-40C 


V(zb)  ’ 


which,  when  substituted  into  (L.52),  yields 


dP 

dz 


dt 

2 

dz 


^CU,z  (1  -<t>)  CUpz 

■f 


V(zb)zb 

V(zb)z? 


V(zb)zb  ) 

(i-<io  cuy ^ 

V(zb)zl 


Differentiating  (L.53)  and  noting  z*  V(z,>),  4>.,  4>,  Ut,  and  Up  are  functions  of  time  only  and  that 
tb  =  Up,  V(zb)  =  AbUp,  and  V(Zb)  =  Ab  we  get 


or 


dP 

dz 


=  -  zC 


4>.fl,  +  4>«t/|  +  (l-4>)  Up  -  ¥Jp 

zbV(zb) 


zC  [4 >♦£/,  +  (l-4>)  Up]  u 

zlV(zb) 

zC  [4>.t7,  +  (l-4>)  Up)  V(zb) 

zbV\zh) 

2 zC  [4>.f J]  +  (l-4>)t/pl 
zlV(zb) 


dp 

dz 


zC 


ZfcV(zb) 


4>A  +  +  (1-WP-WP 


[$>u,  +  (l-$)Up]Up 

Zb 


.Up  +  (1-4 WpW, 
Zb 


zA 


+  2 


PJJ;  +  (i-4»t/;i 


(L.53) 


(L. 54) 


(£,.55) 


48 


or 


ap 

3z 


zC 


hV(zb) 


+  <t >,ut  +  (i-« »0„-$up  + 


z» 


(L.56) 


To  get  (7.  for  substitution  into  (L.56),  we  use  another  form  of  the  solid  phase  momentum  equation 


3k-  dup  go  dP 

IF  +  IF  + 


f. 


p„  dz  (l-e)p. 


(L.57) 


and  noting  that  at  z  =  zb,  dup(zb)/dt  =  and  that  3f/,/3z  =  0  (since  Up  is  a  function  of  time  only),  then 


rr  So  dP  -  _  _ f, _ 

Up  P,  3z  |Z>  (1-e)  p/ 


(L.58) 


Substituting  (L.48)  for  f,  in  (L.58),  we  get 


r,  _  P  W'-Uf  f 
Up  "  d,p,  /,c 


go  3P 

5  S'12- 


Substituting  (L.24)  into  (L.56),  we  get 


(L.  59) 


3/> 

3z 


V,  _  (1^  a 

zbV(z„)  l9*  e  e  '  \ 


|  4,  (*,  ♦  I,  »  -  t/,1  ♦  £  ** 


e  dt 2 


+  t  lir  *  ^  +  (I-*)  - 


2<t>.£/.  ^ 

+  — 


z* 


(L.60) 


7 
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Substituting  (L.58)  into  (L.60),  we  get 


dP 


or 


where 


and 


and 


zC 


80  dz 


z„V(zb) 


t».  -r  + 


a|>*  .  ,  d/np  ri  I  ,  ^  dlnp 

"  7  '  *  'J  e 


**1 
dt2 


<1>X  dlnp  2§.Ug 


dt  ♦  ==*  M  ^  /- 

(L.61) 


a? 
go  az 


<t>2^<,  a?  . 

'  17  7T  1  ‘ 


zC 


'  _  §£  ,y 

zbV(zb)  PP  3z  ‘ 


<t>,  = 


ty.Vp  L§.  dMnp 

4*  -  + 

e  e  dt 2 


(L.  62) 


(L.63) 


<h.E  ^ 

<t>i'  =  b-u,  -  ¥J,  -  -j 


V,  ♦  L  ^  -  U? 


§*L  dlnp  2§.Ut 


e  dt 


+ 


<I>2P 


Zfc 


(Ut-Up)  +  (i Ut-U„T  f, 


(L.(Aa) 


<t>2  = 


(1-e) 


(LMb) 


To  get  £  |  zb  for  elimination  from  (L.62)r  we  evaluate  (L.62)  at  Zj,  and  solve  for  ^|z„.  (L.62) 


evaluated  at  4  is 


dP. 


zbC 


80  dz  IZb  zbV{zb) 


go  dP 
Pp  3z 


r1  ir  I** 


(L.65) 


or 


a? 

dz 


few  {•  -  ;&)’ 


(L.  66) 
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Substituting  (L.66)  back  into  (L.62),  we  get 


or 


or 


where 


t 


z>V(zb) 


QlgoCtyl 

(TT-^Tj 
v  w, 


Pp&oV(.Zb) 


dP 

dz 


zC 

g0zbV(zb) 


1  — — 
v  ***>; 


a? 

dz 


ki(t)z 


*i(0  = 


C 

go*bV(Zb) 


( 


\ 


1- 


4>i 

♦»c 

P/0»> 


/ 


(L. 67) 


(L.68) 


(£..69) 


(L.70) 


kj  depends  upon  the  base  pressure  PB  =  P(Zp)  (one  of  the  values  we  wish  to  solve  for)  through  Vp  (the 
acceleration  of  the  projectile)  and  also  through  (d2lnp)/(dr2). 


Substituting  (L.25), 

^  _  g<ABPB  _  g<ABP, 
'  Mp  Mp 

and  (L.33)  into  (L.70),  we  get 


where 


and 


_ C _  ^'+D-EPra+EPB 

goZbV(Zb)  (\  **c 

V  wJ 


D 


Lfy'Cjjt) 

e 


E 


e 


LA^ 

Vr) 


gpAa 


(L.  71) 


(^•72) 


(L. 73) 
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or 

where 


and 

kl2 

Therefore,  (L.69)  becomes 


K i  —  ^i\Pb  +  ki2 


kn 

ki 


CE^ 

goZhVfa) 


1 


(»,'  +  P) 

goZbVfa) 


k\l  ^ri 


dP 

dz 


(knPB  +  kl2)  z. 


(£.  74) 

(£.  75) 
(£.76) 

(£.77) 

(£.78) 


Integrating  (L.78)  gives  us  the  pressure  distribution  in  the  mixture  0  <  z  <  zb 


P(z)  =  P( 0)  -  -0  ~  z2. 


(£.79) 


We  now  need  the  pressure  distribution  in  the  ullage  region.  For  the  ullage  region,  the  momentum 
equation  for  the  gas  is 

dP 
dz 


0  f  du  du 

go  U  + 


(£.80) 


The  equation  of  continuity  in  the  ullage  region  is 


du 

dz 


1  d£ 

p  dt 


dlnp 
dt  ‘ 


(£.81) 


Integrating  (L.81)  and  noting  (dlnpydt  is  a  function  of  time  only,  we  get 


u(z)  =  - 


dlnp 

dt 


z+m 


(£.82) 
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with  the  boundary  condition, 


«(*,)  =  V  = 


'  +  /<» 


(£•83) 


Therefore 


/«>  -  V,  ♦  *, 


(1.84) 


and 


u(z)  =  - 


d/np  d/np  ,,  .  .  dlnp 

z+Vp+zp  ~dt  = 


dr 


dr 


(£.85) 


Differentiating  (L.85)  gives 


du 

Bt 


...  ,  x  “  lnP 

^  +  (z,~z)  +  i 


dlnp 
”  dt 


(£.86) 


or 


du 

dt 


£  =  K  +  (*,“*>  ^  +  V. 


d2lnp 

dt2 


din  p 

'  df  ‘ 


Substituting  (L.81),  (L.85),  and  (L.87)  into  (L.80),  we  get 


(£.87) 


dP 

dz 


£ 

go 

f 

\ 


Vp+(zp-z) 


d2  lnp 
dt2 


+  V„ 


dlnp 

dt 


v+(z  _z)  ffe' 

p  Kp  '  dtjidt 


or 


ap 

dz 


£ 


£!ep  .  (2_2)  r*sef 

p  {p  Z)  dt2  (  p  ’  {  dt  ) 


(£.88) 


(£.89) 


integrating  (L.89)  from  zb  to  z  and  noting  Zp,  Vp,  and  lnp  are  functions  of  time  only,  we  get 


P(z)  = 


p{h)  ~  7 

6  o 


v  - 


(z,,-z)2  ( d2lnp^ 
dt2 


it,-*)2 


?I)'i 


(£.90) 
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P(z)  =  P(zb) 


v  - 


(Zp-2)2 


where 


A  = 


_  d2lnp  _  ( dlnp 


V 


dt 2  \  dt 


(£•  91) 


(£-92) 


Therefore,  we  get  the  pressure  distribution  in  the  ullage  region  as 

P(z)  =  P(zb)  -  £  (vp(z-zb)  +  [{zp-zbf  -  (zp-zf)  | 


or 


P(z)  =  P(zb)  -  £  [(Vp  +  Zp&)  {z-zb)  -  (z2-z2)  |1 


(£.  93) 


(£.94) 


Defining  the  mean  pressure  Pm  to  be 

f  'p(z)A(z)dz 

P  =  — - 

J  A{z)dz 


and,  since  A(z)  =  AB  is  a  constant,  then 


P  = 


(£.95) 


(£•96) 


Substituting  the  pressure  distribtution  in  the  mixture  region  (L.79)  and  the  pressure  distribution  in  the 
ullage  region  (L.94)  into  (L.96)  to  get  the  mean  pressure,  and  noting  that  P(0)  =  PBr  (the  pressure  at  the 
breech)  we  get 


(knPB+Kn)  2) 

2  2J 


dz 


(z  -ze) 


(£.97) 
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or 


or,  since  L  =  Zp  -  zb, 


r  (k\\P 5  4*  ^12)  a 

=  [P*z*  -  -----  zl  +  P(zt)  (z,-z») 

-  ^((VZ,A>  (Z,-2k)2 

-  |  (z?-z»)  +  Z?A  (Z,-Zt)j 


/*. 


P-  =  p» 


pL 


*11^1 fc2*  ^Uzfr  +  P(zb)L 


^8ozp 


p 

2  f  2 
v*'+ 


(a,  6 zr 


V,+*  AL 
'  3 


) 


(L.98) 


(£•99) 


We  need  to  substitute  for  P^.  To  get  P(zb),  we  evaluate  (L.79)  at  4  and  get 


Substituting  (L.100)  into  (L.99),  we  get 

Pn  = 

Pbt  ~ 


K\ZbP  b 

kxizl 

~~  Pm 

2 

2  * 

p  ^ 

k\\zbP  b 

kqzl 

*  Br 

h 

62 

P 

6z, 

k\\ZbP b 

kazl) 

PL% 

2g0zP 


pL3A 

3g.z,' 


(L.100) 


(L.101) 


Substituting  (L.25)  for  Vp  and  combining  terms,  we  get 

^  *nz>  Lkutl  p£,2Aa  ^ 


P-  =  fV  +  P. 


V 


6z, 


2z,  2zpMp) 


pL?AgP ra  k\2 zl  Lknz\  pL3A 

+  2z,Af,  6zp  2zp  3g0zp' 

Evaluating  (L.94)  at  Zp,  we  get 

P(z,)  =  P(z*)  -  £  [(^+z,A)  (z,-zt)  -  (zp — z^)  |j, 


(L.102) 


(L.103) 
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and  substituting  (L.100)  into  (L.103)  and  noting  L  =  Zp  -  z*,  we  get 


P(zp)  =  PB  =  PBr  - 


+  —  (z2  —  z2)  — 
+  g.  (Z"  h)  2* 


ArjiZfrPt,  fc^z2  pLVp  pLZpA 


80 


go 


Substituting  (L.25)  for  V.,  and  collecting  PB  terms  for  (L.104),  we  get 


0  = 


or 


0  = 


Par + Pa 

l-1- 

pLzpA 

+  £ 

So 

2go 

Par  +  Pa 

-1- 

k\\ZB  ^  pLAB 


kn*l  pLABPtl 


M . 


*n£a 

2 


P^a 

W. 


*ll£b  pL/4BPr< 


M, 


pL2A 
2g»  ‘ 


Subtracting  (L.106)  from  (L.102)  to  eliminate  PPj.  and  get  PB  in  terms  of  Pm,  we  get 


P-  -  P. 


f  knzl  Lk\\tl  pL2AB  4  k\\Z2  pLAB 


6 z, 


p  2zp  2zpMp 


+  1  + 


M, 


p  J 


pL%Pra 

k\%zl 

LkyjZb 

k\2%b 

2 ZpMp 

6zp 

2zp 

+  2 

pLABP ra 

pL2A 

fl—1 

Mp 

2g, 

l 

Noting  the  dependence  of  PB  in  A,  weuse(L.92) 

d2lnp 


A  = 


dt2 


dlnp  Y 

dr  J 


(L.104) 


(L.105) 


(L.106) 


(L.107) 
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or,  substituting  for 


dhnp 
dt 2 


from  (L.33)  and  elminating  Vp  with  (L.25),  we  get 


*  /.n  8<AbPb  ,  8cAbP  , 

A  =  c,(r)  -  + 


vfmp 


mpvf 


(  dlnp 

l  dt 


v 


Substituting  (L.108)  into  (L.107)  and  collecting  PB  terms,  we  get 

knzl  Lknzl  p  L2A, 


6zn 


2zpMp 


+  1  + 


k\\ZB  pLAB 
~2~  + 


AjpL 2 
2  VfMp 


pL2ABPra  p L2  ( 

+  2zpMp  +  2 g, 


r.-fH 

*122* 

pLABP ra 

kXizl 

(7  \ 
Z b  . 

— VL 

^zp  J. 

+  2 

2r 
z,  p 

l3  J 

AjpL2 


1- 


1- 


2 L 
3z 


Pro- 


2L 

3z„ 


-  m) 


pj 


2MpVf 

Substituting  (L.108)  into  (L.106)  and  collecting  PB  terms,  we  get 


PBr  —  Pp 


j  ^  k\\Zb  ^  p LAB 


pL2  f 

+  W.  lc'!,)  " 


dlnp]A 
dt  ) 


A2pL2\ 
2  VrMpj 

AlpL1 
2  VfMp 


knzl  p LAbP r, 


M. 


(L.108) 


(L.109) 


(L.110) 


Equation  (I..  109)  gives  PB  in  terms  of  Pm  (which  is  determined  using  an  equation  of  state)  and  then 
(L.  1 10)  gives  PBr  in  terms  of  PB. 

The  kinetic  energy  (KE)  of  the  gas  and  solid  will  be  required.  The  kinetic  energy  is  given  by 

ga  KE  =  \  //  u2dm  (L111) 
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where  u  is  the  velocity  and  dm  is  the  differential  of  mass;  and  since  pdV  =  dm  and  dV  =  A(z)  dz  =  ABdz, 
then 

1  cXb  1  ctf 

g0  KE  =  -  u2dm  +  2  J,  udm 
or 

g„  KE  =  |  (i-e)  £  u2pppABdz  +  ^  u2pABdz 

+  \  jf'  u2pABdz.  (L.112) 


Substituting  (L.10)  for  Up,  and  (L.14)  for  u  in  the  mixture  region,  and  (L.20)  for  u  in  the  ullage  region,  we 
get 


or 


or 


=  \  <l-e> 

A  H 

PpA$  2 
zb 

I 

**§*¥  I 

z2dz 

► 

*  7  jC 

(v<V 

■*>  7 

j*. 

(L.113) 

g.  M 

ABzb 

6 

[ep  U]  + 

(l-e)  p ,uj] 

oAB  j 

r 

fv,+<b- 

V 

.,)  «*»<>]’  1'. 

;  dr  J  '** 

(L.114) 

KE 

AgZp  2 

=  [eptf;  + 

6 go 

(l-OOpf/p2] 

p AbL 

+  6 go 

w2  + 

3 VJL  ‘""p  +  L* 

'  dr 

(?I! 

..  (L.115) 

It  would  appear  that  the  pressure  used  in  evaluating  the  bunting  rate  should  be  the  mean  pressure  over  the 
region  occupied  by  the  propellant. 
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The  mean  pressure  over  the  mixture  region  is 


P*u 


J  *  P(z)A(z)dz 


f'*  • 

J  A(z)dz 

and  since  A(z)  =  AB  and  using  the  pressure  distribution  in  the  mixture  region  given  by  (L.79),  then 


(L.116) 


^  r 

Ab  J  dz 


(L.  117) 


or 


P*u  -  PBf  ~ 


(knP b  +  kn) 


Zy. 


(L.118) 
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Intentionally  left  blank. 
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APPENDIX  3: 

COMBINED  INFLUENCE  OF  CHAMBRAGE  AND 
PROPELLANT  VELOCITY  LAG 
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Intentionally  left  blank. 
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The  following,  is  work  by  Gough  (Gough,  in  preparation)  with  modifications  and  elaborations  by 
Robbins. 

For  the  effects  of  the  propellant  velocity  lag  and  chambrage  on  the  form  of  the  gradient  equation,  the 
variation  in  area  is  assumed  confined  to  the  chamber  of  the  gun  and  the  propellant  initially  uniformly  fills 
the  chamber. 

The  continuity  equations  for  each  of  the  phases  in  the  mixture  region  are 


9ep  1  BepuA  _ 

i r  *  j  x  ■  '(,) 


5e  l  9  (1— e)  ppUpA 


Pp  Bt  A 


Bz 


=  my 


(*.l) 

(R.2) 


Making  assumptions  analogous  to  the  Lagrange  assumption. 


Bp 

=  0  (gas  density  constant  throughout  the  tube ) 

Be. 

-  0  (porosity  constant  throughout  the  mixture  region) 


(R.  3) 
(R.  4) 


and  noting  the  assumption  that  m)  is  a  function  of  time,  then  (R.1)  becomes, 


or 


3ep  ep  BuA 

IT  +  A  IT  =  m(,)' 


duA  _ 

Bz  ep 


(R.Aa) 


(R.  5) 


Integrating  R.5  and  noting  that  the  right  side  is  a  function  of  time  only  except  for  A  which  is  a  function  of 
z  only,  we  get 


uA 


=  f(t)jAdz  +  k(t )  =  f(t)  V(z)  +  k(t). 


(R.  6) 


Applying  the  boundary  condition  u(0)  =  0  implies  that  k(t)  =  0  since  V(0)  =  0. 
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Also  u(zb)  =  Ug  implies  that 

m 

II 

jsi-s 

w  ^ 
qt 

and,  since  A(zb)  =  AB, 

M(Z) 

U,AbV(z ) 
A(z)V(zby 

Also  (R.2)  becomes  after  using  (R.3)  and  (R.4) 


or 


3e  _  (l-e)Pp  dUpA 
Pp  dt  A  dz 


m 


dUpA  _  w(Q-pp| 

dz  ~  (l-e)pp 


(R.  7) 


(R.  B) 


(K.9) 


(*.10) 


Since  the  right  side  of  (R.10)  is  a  function  of  time  only  except  for  A,  which  is  a  function  of  z  only,  on 
integrating,  we  get 

UpA  =  a  dz  +  ud  =  mv{Z)  +  *,(<).  (/?. id 


Applying  the  boundary  condition,  Up(0)  =  0  implies  that  kt(t)  =  0  since  V(0)  =  0. 
Also  Up(zb)  =  Up  implies  that 


and  since  A(Zb)  =  AB, 


m 


V(zb)  ' 


Up(2) 


UpABV(z) 

A{z)V(zby 


(RA2) 


(R.  13) 


At  the  internal  boundary  defined  by  the  leading  edge  of  the  propellant  bed  we  may  write  the  macroscopic 
balance  of  mass  in  the  following  form 

e  ( U,-Up )  =  Ut+  -  Up  (R.U) 

This  result  includes  the  assumption  that  the  density  jump  across  the  boundary  is  negligible,  as  will 
always  be  the  case  if  the  Mach  number  is  small  compared  with  unity,  and  which  is  in  any  case  consistent 
with  the  present  assumption  R.3. 
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As  with  the  velocity  lag  derivation,  the  ullage  region  is  the  same. 


That  is  the  continuity  equation  given  by 


dp  ,  dp  u 

a  t  dz  ~  ’ 


and  using  (R.3), 


du 

dz 


1  dp 
p  dt 


(R.  15) 


(K.16) 


where  the  right  side  of  (R.16)  is  a  function  of  time  only. 


Therefore, 


“  =  -  p  f  2  +  m 


(R.  17) 


Since  u(z„)  =  Ug+,  then 


y.*  -  -  l  f  2*  +  «»• 


(R.  18) 


and  since  u(Zp)  =  Vp,  then 


v'  '  '  p  f  +  m 


{R.  19) 


Subtracting  (R.18)  from  (R.19)  and  noting 


1  dp 
p  dt 


dfnp 

dt 


we  get 


o’..  -  iir  <*.-*•>  +  v, 


(R. 20) 


Also  solving  (R.19)  for  f3(t)  and  substituting  into  (R.17),  we  get 


u(2)  =  Vp  +  (zp-z) 


(R.  21) 
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Eliminating  Ug+  from  (R.14)  with  (R.20)  gives 


c  (U,-U,)  -  Vr  +  L  ^  -  V, 


or 


u>  -  v> +  i 


V+L  ^ 
P  dt 


-  u. 


(R.  22) 


(R.  23) 


Differentiating  (R.23)  yields 


0*  - 


tL 


1-e 


_  A  f  y  +£  ^?P 

e2  l  '  dr 


-  tf, 


Vp  L  d2lnp  L  dlnp 
+  —  +  -  +  - 


e  dt1 


e  dr 


(K.24) 


We  anticipate  the  need  to  factor  out  the  base  pressure  PB  dependence  of  both  Vp  and  d2lnp/dt2. 


The  projectile  accleration  is  given  by 


V'  *  M. 


8oABP, B  8<AbP r. 


M. 


(R.  25) 


To  get  the  base  pressure  PB  dependence  associated  with  d2lnp/dt2  we  start  with  the  definition  of  the  gas 
density 

m 


P  = 


where  m  =  mass  of  gas  burnt  and  VF  =  free  volume. 


{R. 26) 


Therefore 


and 


m 


V,  =  V„  +  *  (*,-*„>  - 


■e,  A  *  tin  ...  m 

VF  -  ABtp  +  -  ABVp  + 

rp  Vp 


(R.21) 


(R. 28) 
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and 


*  ■  ^  *  l 


where  it  is  understood  that 


may  be  in  the  form 


mm  ,  m 
— »  — .  and  — 
Pp  Pp  Pp 


VW‘  Vm'  -  j 

X — ,  X — ,  and  X — 

P pi  Ppi  Ppi 


where  the  i  index  refers  to  different  propellants. 


{R.  29) 


From  (R.26)  we  get 

and 

and 


inp  =  lnm  -  lnVF 


dlnp  _  rh 
dt  m 


d2lnp 

dt2 


m 

m 


rh2 

m2 


Vp 

vF 


h  + 

v,  vr 


v* 


(fl.30) 


(R.  3D 


(fl.32) 


Substituting  (R.29)  and  then  (R.25),  we  get 


rf2lnp 

dt 


z  =  C,(t) 


8<AbPb  8<AbP i 


where 


Cy(t)  =  m 


m 


vfm„ 


PpVp  j 


MPVP 


n  th2 

+  ~  ~  ~v 
m 


vl 


(R33) 


We  still  require  m,  which  we  get  from 


o  dx 
m  ~  P'S  A 


{R. 34) 
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where  S  is  the  surface  area  of  the  propellant  and  dx/dt  is  the  burning  rate  of  the  propellant  where  dx/dt  is  of 
the  form 


Therefore,  from  (R.34)  we  get 


m 


dx  „« 

—  =  aP  . 
dt  m 


t,  fa  *  d2x 

■  p's  *  +  p-s  ? 


dS  fdxY 

P'  dx  [d*  J 


+  p,S 


£x 

dt2 


and  from  (R.35)  we  get 


~  =  on?:"1 
dr2  dt 


where  dPm/dt  is  determined  numerically  or  from 


/hRT  ^  tRm 
dP m  "w  «V 

~dt  =  V 


—  -  />  ? 


Also,  defining  $  as  the  mass  fraction  of  propellant  burned,  then 


ep  = 


<(>c  -p AbL 
V(zb) 


$JC 

V(z„) 


where 


and 


♦•  =  ♦  ~ 


pAgL 


(1-e)  P,  = 


(1-4»C 

V(zb)  • 


Differentiating  (R.39),  we  get 


=  $  - 


PA 

C 


*  H  ?]• 


(/?.35) 


(R.36) 


(R.  37) 


(R.38) 


(R.39) 


(R.40) 


(*•41) 
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Solving  (R.40)  for  e,  we  get 


e  =  1  - 


(1  -»)C 

P,V(2.) ' 


(RA2) 


Differentiating  (R.42)  and  noting  that 


and  therefore 


V(z„)  =  V(zJ  +  {zb-zpo)AB 


Y(zb)  =  M*  = 


(K.43) 


(rt.44) 


bC  (!-<»  CV(z>)  _  b£  (l-»)  CUjAb 
ppV(zb)  +  p pV\zb)  PPV(2„)  ppV\zb) 


(RA5) 


We  now  have  the  quantiUes  which  arc  required  for  the  solution  of  the  momentum  equation  from  which  we 
get  the  pressure  distrubtion  in  the  mixture  region. 

Consider  the  momentum  equations  for  each  phase 


dAspu  BAepir2 


+  £  «. 


=  -  f,  +  m  up 


(RAG) 


I  f94  (1~£)  +  —  (1~£)  p>“*1  +  (1-e)  go  ~  =  f,  -  rhup  (RA7) 

A  L  dz  j  oz 


where 


/,  =  ^5^  P(u-upffsc 

up 


(K.48) 


fsc  = 


2.5X.  (1-e)  toY 

REN 081  1  (!-£„)  ej  7,°’  ”  ° 


(/?  .49) 
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where 


and 


\  = 


*  £ 
.(>•5 


\2.17 


REN  =  pD,  117,-17,1  /jl 


Adding  (R.46)  and  (R.47)  we  get 


go 


dp 

dz 


-  j  ^  [Aepu  +  A(l-e)p,u,] 
j  £  [Aepu2  +  A(l-e)p,w2]. 


From  (R.8)  we  get 


from  (R.13)  we  get 


from  (R.38)  we  get 


u(z) 


v&m 

A(z)V(zb)' 


Up(z) 


U£AbW) 

A(z)V(zhy 


ep 


<t>^7 

V(z„), 


and  from  (R.40)  we  get 


(1-e)  p. 


(1 HOC 
V(zb) 


which,  when  substituted  into  (R.51),  yields 

$£UtAaV(z)  (1-<| tycu^vjzy 

V\zb)  +  V\zb) 

l  d_  <S>£U]AlV\z)  (1-4>)C77^V2(2)' 

A  dz  V\zb)A{z)  +  V\zb)A(z) 


dP  _  1  d_ 
80  dz  A  dt 


{R.  50) 


(K.51) 


(7?  .52) 
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or 


dP 
80  dt 


CABV(z)  d 
A{z)  dt 

caI  a 


A(z)V\zb)  dz 


Ml  + 

V\zb)  +  V\z„) 


§jj]v\z  )  (Hgjz) 
A(z)  +  A(z) 


(R.  53) 


Performing  the  indicated  differentiation  in  (R.53),  we  get 


dP 
80  dz 


CA,V{z) 

A(z) 

(1 

V\z„) 
CAl 


-  2 


V\zb)  +  V2(zt)  V3(zfc) 


Wp 

v\zb) 


V\zb) 


A{z)V\zb) 


+  (l-W#  ^ 


a) 


(/?.54) 


Noting  that,  since 


V(zb)  =  V(zJ  -»•  /ta  (z*-zj, 


(/?.55) 


and  therefore 

V(zb)  =  -  abup, 

and  since  av(z)/8z  =  A(z),  and  A(z)  and  V(z)  are  function  of  z  only, 


(R  56) 


8z 


V!(z)'|  V\z)  dAiz) 

J  2V(  >  ‘  ^  ’ 


(/?  .57) 


then 


dP 
80  dz 


CAtV(z) 


A(z)V\zb)  L 
+  (1-^)0  p -2 
AbV{z ) 


$.U,  +  2 


<t ).U'AbUp 
V(zb) 


(1-*)^  +  ^  +  (1_^2) 


V{zb) 


V(zb) 


A\z)V(zb) 


(W  +  (l-TO 


8A(z) 

5z 


(K.58) 
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or 


dP 


CABV(z) 


80  dz  A{z)V\zb) 


<1 uOt  +  (1  -0)0,  +  $.U,  -  Wp 


+  2 


ABty»Ug 
V(z„) 

AbV(z) 

A\z)V{zb) 


<U,  -  u,> 


®.u;  +  a -«/,’>  ^  . 


Eliminating  Ut  from  (R.59)  by  using  (R.24),  we  get 


dP 
80  dz 


CABV(z) 

V\z„)A{z) 


<t>«Vp  Lty.  d^p 


+  <J„  1-<1) 


<l>.(l-e) 


-o 


e  dt 2 

v 


!£-vr 


.  L  dlnp 

*•  1  nr 


+  + 


2AB$»Ut 

V(z„) 


(U'-v,) 


AbV(?)  (ATT2  ,  /I  ^(2) 

<W/>  +  °-w')  TT. 


We  need  to  eliminate  (Jp. 


From  (R.2)  we  get,  by  multiplying  by  up, 


de  d  (1-e)  ppUpA 

M'P'  dt  A  dz  ~  Up 


which,  when  substituted  into  (R.47),  gives 


1  dA(\-t)ppup  l  9A(l-e)PpUp  /4  s  dP 

a  — 5 -  +  X  — 5 -  +  (!-£>  *•  ¥ 

.  de  Up  d(l -z)p pUpA 

m  dt  +  A  - fe - 


(R.  59) 


(R.  60) 


(/?61) 


(V?.62) 
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or 


<>-*,  %  *  *=?*  |  ♦  »-«v  £ 


dt  A 

(l-e)p,iij  dA 
A  dz 


„  x  dF 

+  (l-tyo  fa 


=  /,  +  «,  (l-e)pp 


dup  uf(l-e)PpMf  dA 


dz 


dz 


Evaluating  R.63  at  zb  and  noting  dA/dt  =  0, 


and 


dA 

dz 


i,»  =  0  M  2b*Zpo> 


then 


fiP 

(l-e)  8o  fa  +  (l-e)p pVp=f„ 


Using  (R.48)  with  (R.63a),  we  get 


t 1 .  « 


P pdp 


go  dP 

P, 


(R.63) 


{R.63a) 


(R.  64) 
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Substituting  (R.64)  into  (R.60)  gives 


or 


dP 
80  dz 


CAgV(z) 

V\zb)A(z) 


(J).Vp  r  <|>.  d2lnp 

-  +  L —  — 77“ 

.  e  e  dt1 


(t  .  <D-(l-c)l  P (Ut-Up?  , 

f- 


f,  A 

1-<P- 


£1 

p,  a*1'* 


-  ^  fv,  +  i  +  !* 

e2  ^  '  dt  p)  e  dt 

24B<t>.(7. 


V(2,)Az(Z) 


(l  rj2.n  xsr|2v 

3^  [Tr';  »-y.+(,“WP  J 


(/t-6S) 


dP 

9z 


V(z)  CA* 


4(z)  V2(zfc) 


01  — <l>2 


dP 

dz 


«MbV(z)  3A(z) 
V(zb)A\z)  ~dT~ 


(R.  66) 


and 


,  £<!>•  d2Inp  <j>.^p 


t>.  =  +  -T 


e  dr 


(K.67) 


and 


<t>i'  =  <i>.«7,-4>C/„ 


<$>.£ 


p  ^ 


v*  -  V, 


^  dlnp 

j  £  dt 


2AB$,Ut(Ut  -  Up)  .  x  p  (Ut-Up)% 

+  - 7T7-: -  +  ^ - 


V(Z„) 


p PDP 


(R.  68) 


and 


(D,  =  1  -  0  - 


<Ml-e) 


(R.  69) 
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and 


♦,  =  <l>.  t/,1  +  uj. 


(R.  70) 


To  get 


3? 

3z 


for  elimination  from  (R.66),  we  evaluate  (R.66)  at  Zy,  and  solve  for 


With  (R.66)  evaluated  at  Zy,  and  noting 


and  A(zb)  =  Ab,  we  get 


3A 

dz 


'* 


=  0 


or 


dP. 
80  dz 


C 

V(zb) 


(R.  71) 


§R 

dz 


% 


Substituting  (R.72)  back  into  (R.66),  we  get 


or 


V(z)  CAB 
A(z)  v\zb) 


d>i 


PP8oV{ib)  ^1  p^v<^ 


<MjV(z)  8A(z) 
V(z*)A2(z)  3z 


So 


a? 

dz 


V(z)  CM 
Mi)  V\zb) 


4>i 


1- 


p/v*) 


<j>A^(z)  3A(z) 

V(z„)A2(z)  3z 


(R.72) 


(R.73) 


(R.74) 
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<)>!  depends  upon  the  base  pressure  PB  =  P(Zp)  (one  of  the  values  we  wish  to  solve  for)  through  Vp  (the 
accleration  of  the  projectile)  and  also  through  d2lnp/dt2. 

Substituting  (R.33)  into  (R.67),  we  get 


4>. 


A,  M fcC,(0 

+  — : — 


LtykAg  ^ 
eVF  ”  +  e 


(£.75) 


or 


where 


£<j>.  . 

0,  =  V  +  D  +  -f  V, 


(£.  76) 


£4>.Ci(t) 

e 


(£.  77) 


and 


(£.  78) 


Substituting  (R.25)  into  (R.76),  we  get 


<j),  =  (p/  +  D  + 


£<t>. 


B<AbRb  8<AbP  ri 


V  Wp 


A/, 


p  ) 


(£.79) 


or 


!>,  =  «.'  +  £>- 


E^.g0ABPm  E^.g^Ps 


eMp 


eM„ 


(£.80) 


Substituting  «(R.80)  into  (R.74),  we  get 


a? 

dz 


V(z)CAb 

A(z)g0V\zb) 


V  +  D 


<MbV(2)  9A’ 

V(z»)A2(z)  3z 


E**oABFra 


_ 

*2C 

IVVW 


+ 


EU0ABPB 


tMm 


(£.81) 


76 


or 


x  ;  ^  f tf  Et*toA*pa 

a/>  _  cab  _  +  D_  ^T~  ~~^T~  v(z) 

3z  =  g,V2(z>)  1  -  $3  A(z) 

WbV\z)  3/4(2)' 

+  V(z„)A3(z)  3z 


(/?.82) 


Integrating  (R.82)  and  noting  that  A(z)  and  V(z)  are  functions  of  z  and  only  of  z  and  that 


V\z)  3  A 
A3(z)  3z 


3z 


(by  parts)  can  be  written 


V2(z)  + 

2A2(z) 


f  Vfel 
J  A(z) 


dz 


we  get 


P(z) 


P(0) 


C/tB 

g.V2(z») 


<t>,'  +  D 


1 


CAf^  +  flST 

p,*V 


f 

J  A(z) 


<fc 


<MbV2(z)  4A  r  V(z) 
2V(z6>42(z)  +  V(z„)  J  A{z ) 


(R.  83) 
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Noting  P(0)  =  PBr  and  collecting  terms,  we  get 


P(2)  =  Per  + 


cab 

f 

*  ♦  «  -  "X"") 

goV\2„) 

V{zb) 

1  -  .£>  J 

CE^jP, 

v\*Mf 


1  - 


V 


J 


V(z)  4>3C/4bV(z) 

A(z)  2g0V'3(zi))A2(z) 


or 


p,vfr*>  J 

P(z)  =  PBr  +  (alt)  +  alt)  PB)  7,(z)  +  b(t)  Jit) 


(RMa) 


(K.84) 


where 


and 


and 


and 


and 


a1(t) 


CAg 

goV\z„) 


§Ab 

V{zb) 


W  +  D 


1  - 


P/W  ; 


alt) 


cem| 


b(t) 


c<Ml 

2g0V3(z») 


/afr) 


V2(2) 

A2(Z)‘ 


(*.85) 


(K.86) 


(R.%1) 


{RM) 


(R.  89) 


(R.84)  is  the  equation  for  the  pressure  distribution  in  the  mixture  region.  To  complete  the  pressure  distribution 
description,  we  require  the  distribution  in  the  ullage  region  which  should  be  identical  to  the  description  in 
the  velocity  lag  derivation  that  is  in  the  ullage  region.  The  momentum  equation  for  the  gas  is 


dP  p  ( du  du\ 

dz  ~  "  go  U  +  “  fcj* 


(fl.90) 
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The  equation  of  continuity  in  the  ullage  region  is 

8m  1  dp  dlnp 

dz  p  dt  dt 


(R.  91) 


Integrating  (R.91)  and  noting  dlnp/dt  is  a  function  of  time  only  we  get 

dlnp 


u(z) 


dt 


z  +  m 


(/?.  92) 


With  the  boundary  condition  u(zp)  =  V_  we  get 


and  therefore, 


and 


v’  ‘  ~  iir  +  m 


m  -  v,  +  ,,  ^ 


«(z)  =  Vp  +  (z,-z) 


dlnp 
dt  ‘ 


(R.  93) 


(P.94) 


(/?  .95) 


Differentiating  (R.95)  and  noting  zp  =  Vp  gives 

§“  =  *  +  (z  -z)  4!!sp  +  v 
3r  '  }  dr1  '  dt 

Substituting  (R.91),  (R.95),  and  (R.96)  into  (R.90),  we  get 


(R.96) 


dP 

dz 


_P 


V'  +  (z  _z)  +  y 

'  ('  ;  dt 2  dt 


Vp  +  (Z,-2) 


d/np'j  f  dlnp} 
dt  )  dt 


or 


3P 
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P_ 

£0 


,  x  <*2lnP 

+  (z'-z)  ~d^~  '  (2p~z) 


din  pY 
dt  J 


{R.91) 


{R.  98) 


Integrating  (R.98)  from  zb  to  z  and  noting  zp,  and  lnp  are  functions  of  time  only,  we  get 


P(z)  =  P(z6)  - 


v  z  - 


(Zp-zf  d2lnp 


dr2 


(Zp-z)2  r  d/npY 

2  l  dr  Jj 


(P.99) 
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where 


£  V,  -  Jr 

S.  '  2 


rf2Inp 


(/?.100) 


Therefore,  we  get  the  pressure  distribution  in  the  ullage  region  as 


P(z)  =  P(zb) 


-  ±  h 

go  L 


V.  (z-zb)  +  ((z.-zt)  -  (z. 


2  Al 

-2»  2 J 


(/?.102) 


/’(z)  =  P(zb) 


£ 

So 


(V„  +  z, A)  (z-Zj,)  -  (z2-z2) 


(i?.103) 


Defining  the  mean  pressure  Pm  to  be 


/>_  = 


J  />(z)A(z)dz 
J  A(z)dz 


(R.  104) 


and  since#  A(z)dz  =  V(zp)  and  using  (R. 84)  and  (R.  103) 


£.  - 


£*  (Per  +  am  Hz)  +  a&)  PB  Hz)  +  h(l)  Hz))  A(z)dz 


j'f  [p(zb)  -  ft  [(V,  +  Zp  A)  (z-2b)  -  (2 :2-zl )  ~  ])  A(z)dz 


(K.105) 


P_  = 


Par  V(zb)  +  am  Hzb)  +  <*2(0  WO  +  bit)  Uzb) 


+  AgLPfa)  _  Asp  _  J  (23_Z3) 

y(z,)  go  [  V(z,) 


(/?.106) 
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or 


=  P* 


V  (zb)  at(t)  J-iteb)  <h(t)  P$J j(z»)  HO 


where 


and 


V(z„)  V(2p) 


ABLP{zb)  p AB 


V(zp) 


go 


TLA 

V  +  — - 
'  3 


V(zP) 
L 2 


V(z„) 


2  V(zp) 


h{zb)  =  J*  A(z)  /,(z) 


dz 


jm  -  j; 


t/2i 


V\z) 

A(z) 


dz. 


Substituting  in  for  Vp  from  (R.25)  into  (R.107),  we  get 


(R.107) 


(K.108) 


(ft.109) 


P.  =  /*» 


V(z*)  .  ax{tMzb)  a2m{zb)PB  bm(zb) 

+  rr;  ;  +  v  > 


V(zP) 


V(zp) 


AbLP(z„)  pA2BL2PB 


V(zP) 
P AlL2Pra 


V(zp) 


pABL2A 


V(zp) 


2V(zp)Mp  2V{zp)Mp  3  g0V(zpY 


(R.  110) 


To  eliminate  P(Zb),  we  evaluate  (R.84)  at  zb  and  get 

P(h)  -  P„r  +  («i(0  +  ^(0  /*«)  A(zb)  +  b(t)  J2(zb) 


(R.  HD 


and  substitute  into  (R.110)  to  get 


Par 


V(zp) 


ai(t}J'h(zb)P B 
+  V(zp) 


bjtMZb) 

V(zp) 


AsLaJjyj i(zb)  ABLa2(t)J\(zb)PB  ABLb(f)J-^zb) 

+  - rrrr -  +  - nr-: -  +  - mr - 


+  V(zp) 

+ 

V(zp)  4 

V(zp) 

p a\l2pb 

pA2BL2Pra 

pABL2A 

2  V(zp)Mr 

+  2V(z P)MP 

3  g0V(zp) 

Pbr  +  PB  | 

(A.La.mZh) 

l  V(z,) 

pAjL2 

2V(zp)MP 

<h(t)J3{zb) 
f  V{ZP) 

ai(t)h(Zb) 

b{tV<{zb) 

ABLai(t)Ji(zb) 

V(zp)  4 

V(zp)  + 

V(zp ) 

ABLb(t)J1(zb) 

P A2BL2Pra 

pABL2A 

V(zp) 

+  2V(zp)Mp 

3  goV(zpy 

(K.112) 


R.113) 
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We  now  need  PBr  in  terms  of  PB  to  eliminate  PBr  from  (R.1 13)  and  give  us  an  equation  for  PB. 


Evaluating  (R.103)  at  Zp,  we  get 


P(zp)  =  PB 


P{2>) 


(W*) 


L 


or,  since  z2p  -  z2b  =  L  (zp  +  z^, 


Pb 


P<?>) 


&  [v,  + 

8.  V  2  J 


and  substituting  (R.25)  for  Vp,  we  get 


„  „  p LAbPb  ABf)LP ra  PL2A 

PB  ~  P(zb)  ~  — “ —  + 


Af, 


M ’ 


2 go  • 


Substituting  (R.l  1 1)  for  P(Zb)  into  (R.1 16)  and  collecting  terms,  we  get 


0  =  PBr  +  P, 


+  b(tMzb)  + 


-  1  +  ajjyfa)  - 


pLAB 


M, 


+  a, (*)/,&) 


p  / 


Aapf-Pr< 

A/„ 


pL2A 
2g0  ' 


Subtracting  (R.l  17)  from  (R.l  13)  eliminates  PBr  and  gives  PB  in  terms  of  Pm  or 


P_  =  P. 


t  *.\w  /  x  .  pLAB  AMmzb) 

1  -  a^tv^b)  +  -rr-  + 


A/, 


V(zJ 


pAj;.2  ^<V^f*)' 


2V(z,)Af, 


V(zJ 


A'-'  J 


<h  (jW?b)  b(t)J4(zb) 

+  V(zp)  + 


+  AgLa^tMik)  +  ABLb(ty2(zb)  +  pAB2L2ft, 


V(zp) 


V(z.) 


V(Zp) 

pABL3A 
2V(z,W,  "  3g0V(z,) 


-  aWfa)  -  b(t)J2(zb)  -  AbPLP'°  +  P-*A 


A/„ 


2«o 
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(ft.  114) 


(ft. 115) 


(ft.  116) 


(ft.  117) 


(ft. 118) 


Since  A  has  a  dependence  on  PB  through  d2lnp/dt2  using  (R.33)  with  (R.  117),  we  get 


Par  =  Pa 
,  P£ 


1  -  (h{t)Ix{zb)  + 


V 

2  ( 


V 


2  g, 

-  W2&) 

and  using  (R.33)  with  (R.l  18),  we  get 


C,(f)  - 


dlnp 


p LAb  pL2A2B  ^ 

MP  2VfMpJ 

\ 


\  dt  ) 

ABpLP  rt 
Mn 


V  +  g<AaP, 


res 


MpVF  j 


ai(*Vi(zt) 


P.  =  Pn 


1  + 


-  a^y^)  + 

■  3  f 

(. 

\ 


V(zP) 
pLAB 


pABL2  +  a2(t)Ji(zb) 


2  V(zp)MP 


V(zp) 


2  a  2 


pZ/A, 


3/3 


P4JZ, 


P^L 

2&> 


"i 

( 


P  2VFMP 


cm  - 


dlnp 
V  dt  j 


V  + 


M,Vf  , 


3V(z,)VPAfpJ 
1 


2AaL 


W(Zp)j 


aAtMh)  +  b{tMzb)  .  ABLaft)Jfzb) 


V(zP)  V(zp) 

ABLb(t)J2(zb)  p A\L2Pn 


V(zp) 


V(zp) 


2V(zpWP 


-  a.mz,,)  -  b(t)J2(zb)  - 


AppLP f< 

M. 


(7?. 119) 


(/?.120) 


Equation  (R.120)  gives  PB  in  terms  of  Pm  (which  is  determined  using  an  equation  of  state)  and  then  (R.l  19) 
gives  PBr  from  PB. 

The  evaluation  of  Jifo),  J2(zb),  J3(Zb),  and  J4(Zh)  can  be  simplified  by  noting  that  the  variation  in  area 
is  confined  to  the  chamber  and,  therefore. 


and 


J l(Zb)  ~ 


f  a*  *  -  f 

Jo  A(z)  Jo 


for  z 

NJ 

At 

V(z) 

fi?  4- 

c‘i  m 

A(z) 

I4X  T 

A(z) 

dz 


=  +  I 


*  V(zto)  +  AB(z-Zb0) 


dz 


=  JM  +  —  I  V(zJ  (zb-zj  + 


AaiZb-z^y 


.2  A 


(K.121) 


(R.  122) 
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(/?.123) 


,,  %  V2(z),  [V(z*,)  +  Ab  (Zfc-zto)f 

Uh)  m  Wh  = 


Al 


•  r 


V<*>dx] 

A(x)  J 

A(z)dz 

V(jc)dxl 
A(x)  J 

A(z)dz 

f’*{ 

1 

+  A* 

c  (f 


VW 

A(x) 


dx 


ARdz 


AB 


+  J  (z-zj2 


IARdz 


j; 


=  /3(zJ  +  AgJ liz^)  (z0-zj  + 


V(z*>)  (z^-zio)2 


+  g  (Z6  ~  Zjx,)  , 


(K.124) 


« -  r  sg  *  ■  i 


A(z) 


•«w  v^(z) 
A(z) 


dz  + 


I 


'*  (V(zto)  +  Ab  (z-Zi*,)]2 


=  ^(Zi,,,)  + 


tV(zto)  +  A8  (zt-z*)]3  -  V3  (zto) 


3A8 


dz 


(/?.125) 


Equations  (R.122)  -  (R.125)  require  the  evaluation  of  the  integrals  -  J4(Zbo)  only  once. 


The  kinetic  energy  (KE)  of  the  gas  and  solid  will  be  required.  The  kinetic  energy  is  given  by 


go  ke 


u2  dm 


where  u  is  the  velocity  and  dm  is  the  differential  of  mass,  and  since  pdV  =  dm  and  dV  =  A(z)dz,  then 


go  KE  =  \  (1-e)  £  u2ppA(z)dz 

+  |  e  |  u2pA{z)dz  +  ^  f  u2pABdz.  (RA26) 
Z  Jo  Z 


We  get  u  in  the  mixture  region  from  (R.8)  and  up  from  (R.13)  and  u  in  the  ullage  region  from  (R.21),  and 
therefore 


or 


So  KE 


UlAlppV\z) 

V2(zb)A(z) 


+ 


1 

2 


e 


UlAlp„V\z) 

V\zb)A(z) 


+ 


P Ab 
2 


dlnp 

dt 


(Zp~Z) 


(A.  127) 


KE 


(1-e)  UjAfaM  +  afyipJfa) 
2g0V\zb)  +  2g0V\zb) 


pAgL 
6 go 


3Vp  +  3Vp  L 


dlnp  Li  ( dlnpY 
dt  {  dt  )  _ 


(R.  128) 


It  would  appear  that  the  pressure  used  in  evaluating  the  burning  rate  should  be  the  mean  pressure  over 
the  region  occupied  by  the  propellant. 


The  mean  pressure  over  the  mixture  region  is 


and  using  (R.84),  we  get 


P(z)A(z)dz 

I"  A(z)dz 


(R.  129) 


P.  = 

mu 


l  [Par  +  (a,(0  +  a2(t)PB)  Hz)  +  b(t)J2(z)}  A(z)dz 


V(zb) 


or 


p  _  p  .  (fli (O  +  Ths  ,  M.  r,7x 

P—ir  P Br  t //_  \  *3 (Zfr)  "i"  t//_  \  *M\Zfr/* 


V(zb) 


V{zb) 


(R.  130) 


(RA  31) 
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Intentionally  left  blank. 
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APPENDIX  4: 

INPUT  DESCRIPTION  FOR  IBRGA 


87 


Intentionally  left  blank. 
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USER’S  MANUAL  FOR  IBRGA 


IBRGA  relies  on  an  input  data  base  consisting  of  all  numerical  parameters  essential  for  running 
the  code.  All  values  are  in  metric  units.  Below  is  a  compilation  of  a  typical  data  base  showing 
the  name  and  location  of  each  parameter.  The  names  for  the  numerical  values  are  prefixed  with  an 
alphabetical  designator  corresponding  to  the  position  at  which  the  data  are  to  appear,  that  is,  from 
left  to  right  The  data  may  be  separated  by  blanks  or  commas.  The  units  are  shown  to  the  right 
of  each  input 

ABCDEFGHIJK 


record  1 

A.  -  chamber  volume  (cm3) 

B.  -  groove  diameter  (cm) 

C.  -  land  diameter  (cm) 

D.  -  groove/land  ratio  (none) 

(land,  groove,  and  groove/land  ratio  used  to  calculate  the  projectile  bore  area) 

E.  -  twist  (tums/caliber) 

F.  -  projectile  travel  (cm) 

G.  -  gradient  switch  (none) 

(integer  value  designating  the  gradient  equation  to  use 
1  =  Lagrange,  2  =  Chambrage,  3  =  Two-phase,  4  =  RGA) 

H.  -  friction  factor  (none) 

(normally  1  for  granular,  0.01  for  stick,  and  0.1  for  partially  cut  propellant,  only 
used  with  the  two-phase  and  RGA  gradient) 

record  la  (Read  if  and  only  if  gradient  switch  =  2  or  4) 

A.  -  number  of  pairs  of  points  to  describe  chamber  geometry,  integer  I  <=  5  (none) 

B.  -  initial  distance  from  breech  (cm) 

(must  be  0.0) 

C.  -  diameter  at  initial  distance  (cm) 


.  -  Ith  distance  from  breech  (cm) 

(initial  position  of  the  base  of  the  projectile) 

.  -  Ith  diameter  at  Ith  distance  (cm) 

(used  to  calculate  bore  area  overriding  record  1  groove  and  land  diameter 
specifications) 

(note  chamber  geometry  is  used  to  calculate  the  chamber  volume  overriding  record  1 
chamber  volume  description) 


record  2 

A.  -  projectile  mass  (kg) 

B.  -  switch  to  calculate  energy  lost  to  air  resistance, 

integer  (none)  (0  =  no  loss,  1=  loss) 

C.  -  fraction  of  bore  resistance  work  used  to  heat  the  tube  (none)  (>=  0.0,  <=1.0) 

D.  -  gas  pressure  in  front  of  the  projectile  (Pa) 
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record  3 

A.  -  number  of  pairs  of  barrel  resistance  points  (none) 

(integer  J  <=  10) 

B.  -  bore  resistance  (MPa) 

C.  -  travel  (cm) 


.  -  Jth  bore  resistance  (MPa) 

.  -  Jth  travel  (cm) 

record  4 

A.  -  mass  of  recoiling  parts  (kg) 

B.  -  number  of  recoil  point  pairs  (none)  (must  be  an  integer  =  2) 

C.  -  recoil  force  (N)  (force  to  overcome  before  start  of  recoil  -  rod  preload) 

D.  -  time  of  rod  preload  (s)  (must  be  0.0) 

E.  -  recoil  force  (N)  (constant  resistive  force  after  rise  time) 

F.  -  rise  time  (s)  (time  to  go  from  start  of  recoil  to  constant  resistive  recoil  force) 

record  5 

A.  -  free  convective  heat  transfer  coefficient  (W/cm2-K) 

B.  -  chamber  wall  thickness  (cm)  (wall  depth  which  is  heated  uniformly) 

C.  -  heat  capacity  of  chamber  wall  (J/g-K) 

D.  -  initial  temperature  of  chamber  wall  (K) 

E.  -  heat  loss  coefficient  (none)  (usually  1,  but  may  be  set  to  0.0  to  eliminate  heat  loss) 

F.  -  density  of  chamber  wall  (g/cm3) 

record  6 

A.  -  impetus  of  igniter  (J/g) 

B.  -  covolume  of  igniter  (cm3/g) 

C.  -  adiabatic  flame  temperature  of  igniter  (K) 

D.  -  mass  of  igniter  (kg) 

E.  -  ratio  of  specific  heats  of  igniter  (none) 

record  7 

A.  -  number  of  propellants  (none)  (integer  <=  10) 
record  8 

A.  -  impetus  of  propellant  (J/g) 

B.  -  adiabatic  flame  temperature  of  propellant  (K) 

C.  -  covolume  of  propellant  (cm3/g) 

D.  -  mass  of  propellant  (kg) 

E.  -  density  of  propellant  (g/cm3) 

F.  -  ratio  of  specific  heats  of  propellant  (none) 

G.  -  propellant  form  function  indicator  (none) 

(integer,  may  be 

0,  zero  perforated  cylindrical  grain 
1,  one  perforated  cylindrical  grain 
7,  seven  perforated  cylindrical  grain 
15,  nineteen  perforated  hexagonal  grain 
19,  nineteen  perforated  cylindrical  grain) 

H.  -  length  of  propellant  grain  (cm) 

I.  -  diameter  of  perforations  in  the  propellant  grains  (cm)  (value  ignored  if  not  required 

but  must  be  present) 
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J.  -  outside  diameter  of  propellant  grain  (cm) 

(for  the  hexagonal  grain,  it  is  the  distance  between  rounded  comers) 

record  8  repeated  for  each  propellant 

record  9 

A.  -  number  of  burning  rate  triplet  points  (none) 

(integer  J  <=  10) 

B.  -  exponent (none) 

C.  -  coefficient  (cm/s-MPa0) 

D.  -  pressure  (MPa)  (upper  pressure  limit  for  which  the  previous  exponent  and  coefficient 

are  usable) 


.  -  Jth  exponent  (none) 

.  -  Jth  coefficient  (cm/s-MPa”) 

.  -  Jth  pressure  (MPa)  (if  pressure  exceeds  this  limit,  then  this  burning  rate  equation  is 
used) 

record  9  repeated  for  each  propellant 

record  10 

A.  -  integration  time  increment  (ms) 

B.  -  print  increment  (ms) 

C.  -  upper  limit  on  integration  time  to  stop  calculation  (ms) 


* 
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Intentionally  left  blank. 
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APPENDIX  5: 

FORTRAN  LISTING  OF  BRGA 
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Intentionally  left  blank. 
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program  ibrga 

common  nsl  1  ,kpr,fracsl(  1 0),dsdxsl(l  0),surfsl(  10), 

&  nslp(10),tsl(10),pbrch, phase, pmean,bbr(10),abr(10), 

&  deltat,y(20),igrad 
character  outfil*10 
character  bdfile*10 

dimension  br(10),trav(10),tp(10),tr(10),forcp(10),tempp(10),covp< 
&10) 

dimension  chwp(10),ihop(10),gamap(10),nperfs(10),glenp(10),pdp(10 
&),gdiap(10),alpha(10,10),beta(10,10),pres(10,10) 
dimension  a(4),b(4),ak(4),d(20),p(20),z(20),frac(10),surf(10 
&)  ,volp(  1 0),dsdx(  1 0)  jibr(  1 0)  ,ibo(  1 0),tbo(  1 0)  ,d2xdt2(  1 0),tng(  1 0) 
real  lambda j Izp,j2zpj3zpj4zp,jlzbj2zbj3zb,j4zb 
dimension  chdist(5),chdiam(5),bint(4) 
dimension  nsl(  1 0),surfo(  1 0),dsdxn(  1 0) 
c  call  gettim(ihr,imin,isec,ihuns) 

pi=3. 14159 
v.’*ite(*  15) 

15  formatC  input  name  of  data  file  to  be  used  as  input  ’) 
read(*,10)bdfile 

10  format(alO) 

open(unit=2,erp=999,file=bdfile,status=’old’,iostat=ios) 

rewind  2 

write(*,25) 

25  formatf  input  name  of  output  file  ’) 

read(*,10)outfil 

open(unit=6,err=998,file=outfil) 

write(6,16)bdfile 

16  formate  the  input  file  is  ’,alO) 
do  9  i=l,20 

p(i)=0. 

y(i)=o. 

z(i)=0. 

d(i)=0. 

9  continue 

read(2,*,end=20,err=30)cham,grvc, aland, glr,twst,travp,igrad 
&,fsO 

if(igrad.gt.l)go  to  51 
write(6,55) 

55  foimat(lx,’using  Lagrange  pressure  gradient’) 
igrad=l 
go  to  52 

c  define  chambrage  assumes  nctvpts=number  of  points  to  define 

c  chamber  >  or  =  2  <  or  =  5  (?),chdiam(I)  defines  chamber  diameter 

c  at  chdist  (I)  chamber  distance,  chdiam(nchpts)  is  assumed  to  be 

c  the  bore  diameter  and  chdist(i)  is  assumed  to  be  0,  i.s.  at  the 

c  breech.  Assumes  truncated  cones. 

51  if(igrad.eq.3)go  to  401 
if(igrad.eq.4)go  to  434 
write(6,47,err=30) 

47  format(lx, ’Using  chambrage  pressure  gradient’) 
go  to  436 
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434  write(6,437) 

437  foimat(lx, ’using  rga  gradient’) 

go  to  436 

436  read(2  *,end=20,em=30)nchptsI(chdist(I),chdiam(I),I=l,nchpts) 
write(6,53,eir=30)(chdist(I),chdiam(I),I=  1  .nchpts) 

53  format(///,’  chamber  distance  cm  chamber  diameter  cmV(5x,el4 
&.6,5x,el4.6)) 

do  54  1=1  .nchpts 
chdist(I)=0.01  *chdist(I) 

54  chdiam(I)=0.0 1  *chdiam(I) 

c  calculate  chamber  integrals  and  volume 

if(nchpts.gt.5)  write(6,44,err=30) 

44  format(lx,’use  first  5  points’) 
if(nchpts.gt.5)nchpts=5 
bore=chdiam(nchpts) 
if(chdist(l).ne.0.0)write(6,45,err=30) 

45  format(lx,’  #  points  ?  ’) 
chdist(l)=0.0 
pi3=pi/3.0 

bl=0.0 

b2=0.0 

b3=0.0 

b4=0.0 

points=25.0 

56  points=points+points 
step=chdist(nchpts)/points 
zz=0.0 

bint(l)=0.0 

bint(3)=0.0 

bint(4)=0.0 

bvol=0.0 

i2=0.5*chdiam(l) 

k=l 

j=int(points+0.5) 
do  57  1=1  j 
zz=zz+step 

if(k.eq.nchpts-l)go  to  46 
do  58  Il=k,nchpts-1 

if(zz.gt.chdist(I l).and.  zz.ltchdist(U-t-l))go  to  59 

58  continue 
Il=nchpts-1 

59  k=Il 

46  diam=(zz-chdist(k))/(chdist(k+ 1  )-chdist(k)) 
diam=chdiam(k)+diam*(chdiam(k+l)-chdiam(k)) 
rl=0.5*diam 

area=pi*(rl+r2)*(rl+r2)/4. 
bvol=bvol+step*pi3*(rl  *r  1+r  1  *r2+r2*i2) 
bint(  l)=bint(  1  )+step*bvol/area 
bint(3)=bint(3)+step*area*bint(l) 
bint(4)=bint(4)+step,,,bvol*bvol/area 

57  i2=rl 

temp=abs(l  .0-bl/bint(  1)) 

if(abs(  1 .0-b3/bint(3)).gt.temp)temp=abs(l  .0-b3/bint(3)) 
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if(abs(1.0-b4/bint(4)).gt.terap)temp=abs(1.0-b4/bint(4)) 

if(temp.ie.0.001)go  to  41 

bl=bint(l) 

b3=bint(3) 

b4=bint(4) 

go  to  56 

41  cham=bvol*l.e6 
c  write(6,47,err=30)bint(l),bint(3),bint(4) 

c  format(lx,’bint  1  =  ’,614.6,’  bint  3  =  \el4.6,’  bint  4  =  ’,el4. 
c  &6) 

chmlen=chdist(nchpts) 
go  to  52 

401  write(6,402) 

402  format(lx,’using  2  phase  gradient  equation’) 
goto  52 

52  write(6,40,err=30)cham,grve,aland,glr,twst,travp,igrad,fs0 

40  format(lx, ’chamber  volume  cm**3’,el4.6/  groove  diam  cm ’,e  14.67 

&’  land  diam  cm’,el4.6/  groove/Iand  ratio\el4.67’  twist  turns 
&/caliber  *,el4.6 7’  projectile  travel  cm’,el4.6 
&/’  gradient  #  ’43/  friction  factor  ’,el4.6 Ilf) 
cham=cham*l.e-6 
grve=grve*l.e-2 
aland=aland*l.e-2 
travp=travp*l.e-2 

read(2,*,end=20,err^30)prvt,iair,htfr,pgas 

write(6,50,ern=30)prvt,iairdufr,pgas 

50  format(lx, ’projectile  mass  kg’,el4.67’  switch  to  calculate  eneig 

&y  lost  to  air  resistance  J’,i2/  fraction  of  work  against  bore  u 
&sed  to  heat  the  tube’,el4.6/lx,’  gas  pressure  Pa’  ,el4.6) 
read(2,*,end=20,err=30)npts,(br(:),trav(i),i=l4ipts) 
write(6,60,err=30)npts,(br(i),trav(i),i=  1  ,npts) 

60  fomiat(lx, ’number  barrel  resistance  points’, i2/  bore  resistance 

&  MPa  -  travel  cm’/(lx,el4.6,el4.6)) 
write(6,65) 
do  62  i=l,npts 
br(i)=br(i)*l.e6 
trav(i)=trav(i)*  1  .e-2 
62  continue 
65  formatQx) 

read(2,*,end=20,err^30)rcwt,nrp,(rp(i),tr(i),i=  1  jirp) 
write(6,70,ern=30)rcwt,nrp,(rp(i),br(i),i=ljup) 

70  format(lx,’  mass  of  recoiling  parts  kg’,el4.6y’  number  of  recoi 
&1  point  pairs’,i2/  recoil  force  N’,’  recoil  time  sec’/,(lx,el4 
&.6,3x,el4.6)) 
write(6,65) 

read(2,*,end=20,err=30)ho,tshl,cshl,twal,hl,rtiocs 
write(6,75,erp=30)ho,tshl,cshl,twal,hl,rhocs 
75  fomiat(lx,’  free  convective  heat  transfer  coefficient  w/cm**2  K’, 
&el4.6/  chamber  wall  thickness  cm’,el4.6/  heat  capacity  of  st 
&eel  of  chamber  wall  J/g  K’,e  14.67’  initial  temperature  of  chambe 
&r  wall  K’,el4.6/  heat  loss  coefficient’ ,el4.6/  density  of  ch 
&amber  wall  steel  g/cm**3’,el4.6//) 
ho=ho/l.e-4 
tshl=tshl*l.e-2 
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cshl=cshl*l.e+3 
rhocs^rhocs* 1  .e-3/1  .e-6 

read(2,*,end=20,erp=30)forcig,covi,tempi,chwi,gamai 
write(6,85,err=30)forcig,covi, tempi, chwi, gam  ai 

85  format(lx,’  impetus  of  igniter  propellant  J/g’,el4.6 J'  covolume 
&  of  igniter  cm**3/g\el4.6/  adiabatic  flame  temperature  of  igni 
&ter  propellant  K’,el4.6/  initial  mass  of  igniter  kg’,el4.6/  r 
&atio  of  specific  heats  for  igniter’ ,el4.6//) 
forcig=forcig*  1  .e+3 
covi=covi*  1  .e-6/1  .e-3 
read(2,*)nprop 
write(6,98)nprop 

98  formats  there  are  \i2,’  propellants’) 

read(2,*,end=20,err=30)(forcp(i),tempp(i),covp(i),chwp(i), 

&±op(i),gamap(i)mperfs(i),glenp(i),pdp(i),gdiap(i),i=lmprop) 

write(6,95,en^30)(i,forcp(i),tempp(i),covp(i),chwp(i) 

&,rhop(i),gamap(i),nperfs(i),glenp(i),pdp(i),gdiap(i),I=lmprop) 

95  format((’  for  propellant  number’ ,\2J'  impetus  of  propellant  J/g 
&’,el4.6/  adiabatic  temperature  of  propellant  K’,el4.6/  covol 
&ume  of  propellant  cm**3/g’,el4.6/’  initial  mass  of  propellant  kg’ 
&,el4.6 /’  density  of  propellant  g/cm**3’,el4.6/’  ratio  of  specifi 
&c  heats  for  propellant’ ,e  14.6/’  number  of  perforations  of  propell 
&ant’,i2/’  length  of  propellant  grain  cm’,el4.6 /’  diameter  of  per 
&foration  in  propellant  grains  cm’,el4.6 f  outside  diameter  of  pro 
&pellant  grain  cm\el4.6 f)/D 

tmpi=0.0 
do  96  i=l,nprop 
forcp(i)=forcp(i)*  1  .e+3 
covp(i)=covp(i)*  1  .e-6/1  .e-3 
ihop(i)=ihop(i)*  1  .e-3/1  .e-6 

glenp(i)=glenp(i)*0.01 

pdp(i)=pdp(i)*0.01 

gdiap(i>=gdiap(i)*0.01 

tmpi=tmpi+chwp(i) 

kpr=i 

call  prf7 1 0(pdp(i),gdiap(i),glenp(i)jiperfs(i),0. , 

&  fiac(i),volp(i),surf(i),dsdx(i)) 
tng(i)=chwp(i)/rhop(i)/volp(i) 
surfo(i)=surf(i) 
write(6,408)i,tng(i) 

408  format(’  for  propellant  \i2,’  the  total  number  of  grains’ 

&  ,’  is  \el4.6) 

96  continue 
tmpi=tmpi+chwi 
do  97  j=l,nprop 

read(2,*,end=20,err=30)nbr(j),(alpha(j,i),beta(j,i),pres(j,i), 

&i=l,nbr(j)) 

write(6, 1 1 0,err=30)nbr(j),(alpha(j  ,i)  ,beta(j  ,i),pres(j  ,i), 

&i=l,nbr(j)) 

110  format(lx, ’number  of  burning  rate  points ’,i2/3x,’  exponent’,8x,’ 

&  coefficient’.lOx,’  pressure’/5x,’-’,15x,’cm/sec-MPa**ai’,10x,’MP 
&a  V(  1  x,e  14.6^x,e  14.6, 1 5x,el4.6)) 
do  112  i=l,nbr(j) 
beta(j»i)=beta(j,i)*  1  .e-2 
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presG,i)=pres(j,i)*  1  e6 

1 12  continue 

97  continue 
writc(6,65) 

rcad(2,*,end=20,ern=30)deltat,deltap,tstop 

wiitc(6,120,cn^30)deltat,deltap,tstop 

120  format(lx,’time  increment  msec’,el4.6,’  print  increment  msec’,el4 
&.6/lx,’time  to  stop  calculation  msec  \el4.6) 
write(*,130) 
deltat=deltat*0.001 
deltap=deltap*0.00 1 
tstop=tstop*.001 

130  format(lx,’end  input  data  --  LB.  calculation  start’) 
if(igrad.eq.2.or.igrad.eq.4)go  to  131 
borc=(glr*grve*grve+aland*aland)/(glr+ 1 .) 
bore=sqrt(bore) 

131  areab=pi*bore*bore/4. 

lambda= 1 ./(( 1 3.244.  *log  1 0(1 00.  *bore))**2) 

iplot=0 

pltdt=deltat 

pltt=0. 

pmaxm=0.0 

pmaxbMJ.O 

pmaxba=0.0 

tpmaxm=0.0 

tpmxbr=0.0 

tpmxba=0.0 

tpmax=0.0 

a(l)=0.5 

a(2)=l.-sqrt(2.)/2. 

a(3)=l.+sqrt(2.)/2. 

a(4)=l./6. 

b(l)=2. 

b(2)=l. 

b(3)=l. 

b(4)=2. 

ak(l)=0.5 

ak(2)=a(2) 

ak(3)=a(3) 

ak(4)=0.5 

vp0=0.0 

tib=0.0 

tcw=0.0 

if(igrad.eq.3)chmlen=cham/areab 

zb=chmlen 

zp=chmlen 

grlen=0. 

grdiam=0. 

egama=0. 

do  5  i=l,nprop 

grlen=grien+chwp(i)*glenp(i) 

grdiam=grdiam+chwp(i)*gdiap(i) 

ibo(i)=0 

egama=egama+chwp(i)*gamap(i) 
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nsl(i)=0 

5  vpO=chwpO)/AopO)+vpO 
volgi=cham-vpO-chwi*covi 
giien=giien/(tmpi-chwi) 
grdiam=grdiam/(tmpi-chwi) 
egama=(egama+chwi*gamai)Ampi 
ism=0 
odlnr=0. 
vfO=cham-vpO 
epsO=l.-vpO/cham 
eps=epsO 
gasden=chwi/vlO 
prden=tmpi/vpO 
ug=0. 
up=0. 

pmean=forcig*chwi/volgi 

pbase=pmean 

pbrch=pmean 

opbase=pmean 

volg=volgi 

volgi=volgi+vpO 

wallt=twal 

tgas=tempi 

told=0. 

tgaso=tgas 

dtgaso=0. 

covl=covi 

t=0. 

ptime=0.0 

ibrp=10 

z(3)=l. 

nde=ibrp+nprop 

wri  te(6, 1 32)areab,pmean,vpO,volgi 

132  foimat(lx,’area  bore  mA2  \el 6.6/  pressure  from  ign  pa’,el6.6y 
volume  of  unbumt  prop  mA3  ’,el6.6 J 
&*  init  cham  vol-cov  ign  mA3  \el6.6) 
write(6,6) 

6  format(lx,’  time  acc  vel  dis  mpress 
&  pbase  pbrch  ’) 
iswl=0 

19  continue 
do  11  J=l,4 

c  FIND  BARREL  RESISTANCE 

do  201  k=2,npts 
if(y(2)+y(7).ge.trav(k))go  to  201 
go  to  203 

201  continue 
k=npts 

203  rcsp=(trav(k)-y(2)-y(7))/(trav(k)-trav(k-l)) 
resp=br(k)-resp*(br(k)-br(k- 1 )) 

c  FIND  MASS  FRACTION  BURNED 

do  211  k=l,nprop 
kpr=k 

if(ibo(k).eq.  l)goto21 1 
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nsll=0 

call  prf7 10(pdp(k),gdiap(k),glenp(k))nperfs(k),y(ibrp+k) 
&,frac(k),volp(k),surf(k),dsdx(k)) 
nsl(k)=nsll 

if(nsl(k).eq.O)goto  212 
if(nslp(k).eq.l)go  to  212 
write(6,2213)k 

2213  format(’  propellant’, i2,’  has  slivered’) 
nslp(k)=l 
tsl(k)=y(3) 
ism=l 

212  continue 

if(frac(k).lt..9999)  go  to  211 

frac(k)=l. 

tbo(k)=y(3) 

ibo(k)=l 

ism=l 

write(6,456)k 

456  formate  propellant’ ,i2,’  has  burned  out’) 

211  continue 

c  ENERGY  LOSS  TO  PROJECTILE  TRANSLATION 
elpt=prwt*y(  1  )*y(  1  )/2. 
eptdot=prwt*y  ( 1 )  *z(  1 ) 

c  ENERGY  LOSS  DUE  TO  PROJECTILE  ROTATION 
elpr=pi*pi*prwt*y(l)*y(  1  )/4.*twst*twst 
eprdot=pi*pi*prwt*y(l)*z(l)/2.*twst*twst 
c  ENERGY  LOSS  DUE  TO  GAS  AND  PROPELLANT  MOTION 
if(igrad.eq.l)go  to  214 
if(igrad.eq.3)goto  217 
if(igrad.eq.4)go  to  438 
pt=y(2)+y(7) 
vzp=bvol+areab*pt 

j4zp=bint(4)+((bvol+areab*pt)**3-bvol**3)/3./areab/areab 

elgpm=tmpi*y(l)*y(l)*areab*areab*j4zp/2./vzp/vzp/vzp 
go  to  216 

438  pb=y(7)+y(10) 

vzb=bvol+areab*pb 

j4zb=bint(4)+(vzb*  *  3-bvol*  *3)/3  ./areab/areab 
elgpm=(  1  .-eps)*up*up*areab**2*piden*j4zb+ 

&  eps*ug*ug*areab**2*gasden*j4zb 
elgpm=elgpm/2./vzb/vzb+gasden*areab*'ullen/6.* 

&  (3.*y(l)*y(l)+3.*y(l)*ullen*dlnrho+ullen**2*dlniho**2) 
c  approximate  epdot 

epdot=tmpi*y(l)*z(l)/3. 
go  to  216 

214  elgpm=tmpi*y(l)*y(l)/6. 
go  to  216 

217  elgpm=areab*zb/6.*(eps*gasden*ug*ug+(l.-eps)*prden*up*up) 
elgpm=elgpm+gasden*areab*ullen/6  *(3.*y(  1  )*y(  1 )+ 

&  3.*y(l)*ullen*dlnrho+ullen**2*dlnrho**2) 
c  approximate  epdot 

epdot=tmpi*y(  l)*z(  l)/3. 

c  ENERGY  LOSS  FROM  BORE  RESISTANCE 
216  elbn=y(4) 
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z(4)=areab*resp*y(l) 

ebrdot=z(4) 

c  ENERGY  LOSS  DUE  TO  RECOIL 
elrc=rcwt*y(6)*y(6)/2. 
erdot=rcwt*y(6)*z(6) 

c  ENERGY  LOSS  DUE  TO  HEAT  LOSS 

areaw=cham/areab*pi*bore+2.*areab+pi*bore*(y(2)+y(7)) 

avden=0.0 

avc=0.0 

avcp=0.0 

zl8=0 

zl9=0 

do  213  k=l,nprop 

zl  8=forcp(k)*gamap(k)*chwp(k)*frac(k)/(gamap(k)- 1  ,)/tempp(k)+zl  8 

zl  9=chwp(k)*frac(k)+z  19 

avden=avden+chwp(k)*frac(k) 

213  continue 

avcp=<zl8+forcig*gamai*chwi/(gamai-L)/tempi)/(zl9+chwi) 

avden=(avden+chwi)/(volg+cov  1 ) 

awel=.5,‘y(l) 

htns=lambda*avcp*avden*awel+ho 

z(5)=areaw’,‘htns*(tgas-wallt)*hl 

elht=y(5) 

ehdot=z(5) 

wallt=(elht+htfr,'‘elbr)/cshl/rhocs/areawAshl+twal 
c  write(6,*)lambda,avcp,avden,awel,ho,areaw,htns,tgas,wallt,hl,z(5) 
c  &,elht 

c  ENERGY  LOSS  DUE  TO  AIR  RESISTANCE 

ain^iair 

z(8)=y(l)*pgas*air 
elar=areab*y(8) 
eddot=z(8)*areab 
c  RECOIL 

z(6)=0.0 

if(pbrch.le.rp(l)/areab)go  to  221 
rfor=rp(2) 

if(y(3)-ti0.ge.tr(2))go  to  222 

rfor=(tr(2)-(y(3)-tiO))/(tr(2)-tr(l)) 

rfor=rp(2)-rfor*(ip(2)-ip(l)) 

222  z(6)=areab/rcwt*(pbrch-rfor/areab-resp) 
if(y(6).lt.0.0)y(6)=0.0 

z(7)=y(6) 
goto  223 
221  tri>=y(3) 

223  continue 

c  CALCULATE  GAS  TEMPERATURE 
eprop=0.0 
rprop=0.0 
dmfogt=0.0 
dmfog=0.0 
do  231  k=l,nprop 

eprop=eprop+forcp(k)*chwp(k)*frac(k)/(gamap(k)- 1 .) 
iprop=rprop+forcp(k)*chwp(k)*frac(k)/(gamap(k)- 1  .)Acmpp(k) 
dmfogt=dmfogt+forq>(k)*iiiop(k)*tng(k)*surf(k),,,z(ibrp+k)/ 
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&  ((gamap(k)-L)*tempp(k)) 

dmfog=dmfog+forcp(k)*rhop(k)*tng(k)*surf(k)*z(ibrpfk)/ 

&  (gamap(k)-l.) 

231  continue 

tenerg=elpt+elpr+elgpm+elbr+elrc+elht+elar 
tgas=(eprop+forcig*chwi/(gamai- 1  .)-elpt-elpr-elgpm-elbr-elrc-elht 
&-elar)/(tprop+forcig*chwi/(gamai- 1  .)/tempi) 
tedot=epdot+eprdot+eddot+ebrdot+erdot+ehdot+eptdot 
dtgas=(dmfog-tedot-tgas*dmfogt)/(rpropfforcig*chwi/ 

&  (gamai-l.)/tempi) 
c  FIND  FREE  VOLUME 

v  1=0.0 
cov  1=0.0 
do  241  k=l,nprop 
v  l=chwp(k)*(  1  .-frac(k))/ihop(k)+v  1 
cov  1  =cov  1 +chwp(k)  *covp(k)*frac(k) 

241  continue 

volg=volgi+areab*(y(2)+y(7))-v  1  -cov  1 
c  CALCULATE  MEAN  PRESSURE 

r  1=0.0 

do  251  k=l,nprop 

r  1 =r  1 +forcp(k)*chwp(k)*frac(k)/tempp(k) 

251  continue 

pmean=tgas/volg*(rl+forcig*chwi/tempi) 

259  resp=resp+pgas*air 
if(igrad.eq.l)go  to  252 
if(igrad.eq.2)goto  403 
if(igrad.eq.3)go  to  404 
if(igrad.eq.4)go  to  441 

403  if(iswl.ne.0)go  to  253 
pbase=pmean 
pbrch=pmean 
if(pbase.gt.resp+l.)iswl=l 
go  to  257 

c  USE  CHAMBRAGE  PRESSURE  GRADIENT  EQUATION 
253  j  1  zp=bint(  1  )+(bvol*pt+areab/2.  *pt*pt)/areab 
j2zp=(bvol+areab*pt)**2/areab/areab 
j3zp=bint(3)+areab*bint(l)*pt+bvol*pt*pt/2.+arcab*pt'l'pt*pt/6. 
a2t=-tmpi*areab*areao/prwt/vzpA'zp 
alf=l.-a2t*jlzp 

alt=tmpi*areab*(areab*y(l)*y(l)/vzp+areab*resp/prwt)/vzp/vzp 
bt=-tmpi*y(  1  )*y(l  )"'arcab*areab/2./vzp/vzpA,zp 
bata=-alt*j  lzp-bt*j2zp 
gamma=alf+a2t*j3zp/vzp 
delta=bata+a  1  t*j  3zp/vzp+bt*j4zp/vzp 
c  calculate  base  pressure 

pbase=(pmean-delta)/gamma 
c  calculate  breech  pressure 

pbrch=alf*pbase+bata 
go  to  254 

c  USE  2  PHASE  GRADIENT  EQUATION 

404  IF(ISW1.NE.O)GOTO  407 
pbase=pmean 
pbrch=pmean 
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if(pbase.gtresp+ 1  )isw  1 = 1 
go  to  257 

407  if(iswl.eq.2)go  to  411 

vzp=cham+arcab*(y(2)+y(7)) 

vzb=cham+areab*(y(10)+y(7)) 

phi=0. 

phidot=0. 

dmorho=0. 

dmcov=0. 

dmromw=0. 

rmomw=0. 

vfree=vzp-vl 

do  405  k=l,nprop 

rmomw=nnomw+chwp(k)*frac(k)*forcpOt)Aempp(k; 

phi=chwp(k)*frac(k)+phi 

if(ibo(k).eq.l)go  to  405 

dmorho=dmoiho+tng(k)*surf(k)*z(ibrpfk) 

phidot=rhop(k)*tng(k)*surf(k)*z(ibrp+k)+phidot 

dmcov=ihop(k)*tng(k)*surf(k)*z(ibrp+k)*covp(k)+dmcov 

dmromw=dmromw+iiiop(k)*tng(k)*surf(k)’,'z(ibrp+k)* 

&  forcp(k)/tempp(k) 

405  continue 

nnomw=rmomw+chwi*forcig/tempi 

gasmas=phi+chwi 

gasden=gasmas/vfree 

phi=(phi+chwi)Ampi 

if  (phi.gt.0.999)  then 

iswl=2 

rbm=pbase/pmean 
rbrm=pbrch/pmean 
if(phi.ge.l.)go  to  411 
endif 

dradt=phidot 
phidot=phidot/tmpi 
vdotov=(dmoiho+areab*y(l))/vfree 
dlnrho=dmdt/gasmas-vdotov 
dvoldt=dmorho+areab*y(  l)-dmcov 

c  GET  TIME  DERIVATIVE  OF  MEAN  PRESSURE 

dpmdt=(dmromw*tgas-pmean*dvoldt+dtgas*nnomw)Arolg 

volprp=0. 

effdia=0. 

dmdmdt=0. 

dmdmop=0. 

avelen=0. 

avedia=0. 

do  406  k=  laptop 

if(ibo(k).eq.l)go  to  406 

volprp=volprp+(l  .-frac(k))*chwp(k)/rhop(k) 

dmdmdt=dmdmdt+rhop(k)*mg(k)*dsdxOO*z(ibrp+k)*z(ibrp+k) 

dmdmdt=dmdmdt+rt»p(k)*tng(k)*surfi(k)1',d2xdt2(k) 

dmdmor^dmdmor+(dsdx0c)*z(ibip+k)**2+surf(k)*d2xdt2(k))*tng(k) 

effdia=effdia+6.*volp(k)/surf(k)*(l.-frac(k))*chwp(k) 

406  continue 

c  1  t=dmdmdt/gasmas-dmdmor/vfree+vdotov*  *2-(dmdt/gasmas)*  *2 
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d21nr=clt-areab**2*pbase/vfree/prwt 

d21nj^d21nr+areab*areab*resp/vfrce/prwt 

zp=chmlen+y(2)+y(7) 

zb=chmlen+y  ( 1 0)+y  (7) 

ullen=zp-zb 

cnow=tmpi-gasmas 

vp=y(i) 

effdia=effdia/cnow 

prden=cnow/volpip 

up=y(9) 

phistn=phi-gasden*areab*ullen/tinpi 

ulldot=vp-up 

dphist=phidot-gasden*areab/tmpi*(ulldot+ullen*dlniho) 
eps= 1  .-(l.-phi)*tmpi/prden/vzb 
epsdot=phidot*tmpi/prden/vzb+(l.-phi)*tmpi*up*areab/ 

&  piden/vzb/vzb 
ug=up+(vp+ullen*dlnrtio-up)/eps 
alam=(l  .5*grlen/grdiam)**.666666667 
alam=(0.5+grlen/grdiam)/alam 
alam=alam**2.17 
c  VIS  kg/sAn 
vis=. 00007 

ren=gasdcn/vis*effdia*abs(ug-up) 

if(ren.lt.l.)ren=l. 

fsrg=2.5*alam/ren**.08 1  *((1  .-eps)/(l.-eps0)*eps0/eps)**.45 
fsc=fsrg*fsO 

phi2=  l.-phi-phistr*(l  .-eps)/eps 
phi  1  p=dphist*ug-phidot*up-phistr*epsdot/eps/eps 
&  *(vp+ullen*dlnrho-up)+pbistr,,,ulldot*dlnriK)/eps 
&  +2.*phistr*ug/zb*(ug-up) 

phi  lp=plulp+phi2*gasden/effdia/prden*(ug-up)**2*fsc 
ak2=  1 ./( 1 .  -phi2*tmpi/prden/vzb) 
phil=phi  1  p+phistr*z(  l)/eps+ullen*phistr*d21nr/eps 
c  ACCELERATION  OF  FORWARD  BOUNDARY  OF  PROPELLANT  BED 
z(9)=gasden*(ug-up)**2*fs{^rden/effdia+tmpi*phil*ak2 
&/vzb/prden 
z(10)=y(9) 

e=phistr/eps*(  1  .-uUen*areab/vfree)*areab/prwt 
dd=ullen*phistr*c  1  t/eps 
akl  l=tmpi*e*ak2/zb/vzb 
akl2=tmpi*ak2*(philp+dd)/zb/vzb-akl  l*resp 
pbase=pmean-akl2*zb*zb/2.+gasden*uUen*resp*areab/prwt 
pbase=pbase+ak  1 2*zb*zb*(zb/3  .+uiIen)/2./zp 
pbase=pbase-gasden*ullen**2*areab*resp/2./zp/prwt 
pbase=pbase-gasden*ullen**2/2.*(l.-2.*ullen/3./zp)* 

&  (clt-dlmho**2) 

pbase=pbase-areab’,'*2*gasden*ullen**2* 
&(l.-2.*uUeiV37zp)*resp/prwt/vfree/2. 
deno=-akl  l*zb**3/6./zp-ulIen*akl  l*zb*zb/2./zp 
deno=deno+gasden*ullen*areab/prwt-areab**2*gasden*ullen**2 
&*(1  .-2.*uUen/3./zp)/2./vfree/prwt 

deno=deno-gasden*ullen**2*areab/2./zp/prwt+l.+akll*zb*zb/2. 

pbase=pbase/deno 

if(ism.eq.0)goto453 
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ifi(ism.eq.l)goto451 

gOto452 

451  ism=2 

tss=sqrt(egama*rmomw/gasmas*tgas) 

write(6,*)tss 

tss=ullen/(ullen*odlnr+tss) 

tso=y(3) 

write(6,*)tss,tso 

452  coefbp=(tss+tso-y(3)-deltat)/tss 
if(coefbp.gt.  1  .)coefbp=  1 . 
if(coefbp.le.O.)then 
coefbp=0. 

ism=0 

endif 

pbase=coefbp*opbase+(l  ,-coefbp)*pbase 
write(6,*)coefbp,opbase,pbase,ism 

453  odlnr=d]nriio 
opbase=pbase 

pbrch=pbase*(l.+akl  l*zb*zb/2.+gasden*ullen*areab/prwt 
&  -areab**2*ga?den*ullen**2/2./vfree/prwt) 
pbrch=pbrch+ak;2*zb*zb/2.-gasden*ullen*areab*resp/prwt 
pbrch=pbrch+gasden*ullen**2/2  *(clt-dlniiio**2) 
pbrch=pbrch+areab**2*gasden*ullen**2*resp/2./vfree/prwt 

go  to  254 

C  USING  RGA  GRADIENT 

441  if(iswl.ne.O)go  to  444 
pbase=pmean 
pbrch=pmean 
if(pbase.gLresp+l.)»swl=l 
go  to  257 

444  if(iswl.eq.2)go  to  411 

vzp=cham+areab*(y(2)+y(7)) 
vzb=cham+areab*(y(10)+y(7)) 
j  1  zb=bint(  l)+(bvol*pb+areab/2.*pb*pb)/areab 
j2zb=(bvol+areab*pb)  *  *2/areab/areab 

j3zb=bint(3)+areab*bint(l)*pb+bvol*pb*pb/2.+areab/6.*pb**3 

jrfil=0. 

phidot=0. 

dmorho=0. 

dmcov=0. 

dmromw=0. 

nnomw=0. 

vftee=vzp-vl 

do  442  k=l,nprop 

imomw=rmomw+chwp(k)*frac(k)*forcp(k)/tempp(k) 

phi=chwp(k)*frac(k)+phi 

if(ito(k).eq.l)go  to  442 

dmorho=dmorbo+tng(k)*surf(k)*z(ibrp+k) 

phidot=ihop(k)*tng(k)*suif(k)*z(ibrpH-k)+phidot 

dmcov=rtiop(k)*tng(k)*surf(k)*z(ibip+k)*covp(k)+dincov 

dmromw=dmromw+rhop(k)*tng(k)*surf(k)*z(ibrp+k)* 

&  forcp(k)/tempp(k) 

442  continue 

rmomw=rmomw+chwi*forcig/tempi 
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gasmas=phi+chwi 
gasden=gasmas/vfree 
phi=(phi+chwi)/tmpi 
if  (phi.gt.0.99)  then 
iswl=2 

rbm=pbase/pmean 
rbrm=pbrch/pmean 
if(phi.ge.l.)go  o  411 
endif 

dmdt=phidot 
phidot=phidot/tmpi 
vdotov=(dmoriio+areab*y(l))/vfree 
dlmho=dmdt/gasmas-vdotov 
dvoldt=dmorho+areab*y(  1  )-dmcov 
c  get  time  derivative  of  mean  pressure 

dpmdt=(dmromw*tgas-pmean*dvoldt+dtgas*rmomw)/volg 

volprp=0. 

effdia=0. 

dmdmdt=0. 

dmdmor=0. 

avelen=0. 

avedia=0. 

do  443  k=l,nprop 

if(ibo(k).eq.l)go  to  443 

volprp=volprp+(l  .-firac(k))*chwp(k)/rhop(k) 

dmdmdt=dmdmdt+rhop(k)*mg(k)*dsdxOc)*z(ibrp+k)*z(ibrp+k) 

dmdmdt=dmdmdt+rhop(k)*mg(k)  *surf(k)*d2xdt2(k) 

dmdmoi^dmdmor+(dsdx(k)*z(ibrp+k)'*,*2+surf(k)*d2xdt2(k))*tng(k) 

effdia=effdia+6.*volp(k)/surf(k),,'(l.-fi3c(k))*chwp(k) 

443  continue 

c  1  t=dmdmdt/gasmas-dmdmor/vfree+vdotov*  *2-(dmdt/gasmas)*  *2 

d21nr=c  1  t-arcab  *  *2*pbase/vfree/prwt 

d21nr=d21nr+areab*areab*resp/vfree/prwt 

zp=chmlen+y(2)+y(7) 

zb=chmlen+y(  1 0)+y  (7) 

ullen=zp-zb 

cnow=tmpi-gasmas 

vp=y(l) 

effdia=effdia/cnow 

prden=cnow/volprp 

up=y(9) 

phistp=phi-gasden*areab*ullen/tmpi 

ulldot=vp-up 

dphist=phidot-gasden*areab/tmpi*(ulldot+ullen*dlnrho) 
eps=  1 .-( 1  .-phi)*tmpi/prden/vzb 
epsdot=phidot*tmpi43rden/vzb+(l.-phi)*tmpi*up*areab/ 

&  prden/vzb/vzb 
ug=upf(vp+ullen*dlnrho-up)/eps 
alam=(1.5*grlen/grdiam)**.666666667 
alam=(0.5+grlen/grdiam)/alam 
alam=alam**2.17 
c  VIS  kg/s/m 
vis=.00007 

ren=gasden/vis*effdia*abs(ug-up) 
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if(ren.lt.l.)ren=l. 

fsrg=2.5*alam/ren**.08 1  *((1  .-eps)/(l  .-epsO)*cpsO/eps)**.45 
fsc=fsrg*fsO 

phi2=l.-phi-phistr*(l  ,-eps)/eps 
philp=dphist*ug-phidot*up-phistr*epsdot/eps/eps 
&  *(vp+uhen*dlnrho-up)+phistr*ulluoi*dlnrtio/eps 
&  +2.  *areab*phistr*ug/vzb*(ug-up) 
philp=philp+phi2*gasden/effdia/prden*(ug-up)**2*fsc 
ak2=  1 ./( 1 .  -phi2*tmpi/prdeiVvzb) 
phi  l=phi  lp+phistr*z(  l)/eps+ullen*phistr*d21nr/eps 
c  ACCELERATION  OF  FORWARD  BOUNDARY  OF  PROPELLANT  BED 
z(9)=gasden*(ug-up)**2*fsc/prden/effdia+tmpi,!iphil*ak2 
&/vzb/prden 
z(10)=y(9) 

phi3=phistr*ug*ug+(l  .-phi)*up*up 
e=  1  .-ullen*areab/vfree 
dd=ullen*phistr*clt/eps 

alt=tmpi*areab/vzb/vzb*(phi3*areab/vzb-(philp+dd-e 
&  *phistr*areab*resp/eps/prwt)*ak2) 
a2t=(-tmpi*e*phistr*areab**2/vzb/vzb/eps/prwt)*ak2 
bt=-tmpi*phi3*areab**2/2./vzb/vzb/vzb 
pbase=pmean-gasden*ullen**2/2.*(c  1  t-dlnrho**2 
$  +areab**2*resp/prwt/vfree)*(l  .-2.*areab*ullen/3./vzp) 
pbase=pbase-al  t*j3zb/vzp-bt*j4zb/vzp-areab*ullen*a  l  t*j  1  zb/vzp 
pbase=pbase-areab*bt*uUen*j2zb/vzp-gasden*aieab**2*ullen**2 
&*rcsp/2./vzp/prwt+alt*jlzb4-bt*j2ztH-areab*gasden*ullen*resp/prwt 
deno=l.+arcab*ullen*a2t*jizbAzp-gasdcn*arcab**2*ullen**2 
&  /2./vzp/prwt+a2t*j3zb/vzp-a2t*j  1  zb+gasden*ullen*areab/prwt 
deno=deno-gasden*ullen**2*areab**2/2./vfrce/prwt 
&  +gasden*areab**3*ullen**3/3./vzp/vfree/prwt 
pbase=pbase/deno 

pbrch=pbase*(l  .-a2t*j  lzbf  gasden*ullen*areab/prwt 
&  -gasden*ullen*  *2*areab*  *2/2  .A  free/prwt) 
&+gasden*uUen**2/2.*(clt-dlniho**2+areab**2*resp/prwt/vfree) 

&  -alt*j  lzb-bt*j2zb-areab*gasden*ullen*resp/prwt 
go  to  254 

411  pbase=rbm*pmean 
pbrch=rbrm*pmean 
go  to  254 

c  USE  LAGRANGE  PRESSURE  GRADIENT  EQUATION 
252  if(iswl.ne.O)go  to  256 

if(pmean.lLresp)resp=pmean 
c  CALCULATE  BASE  PRESSURE 

256  pbase=(pmean+tmpi*resp/3./prwt)/fl  .+tmpi/3./prwt) 
if(pbase.gt.resp+l.)iswl=l 

c  CALCULATE  BREECH  PESSURE 
pbrch-pbase+tmpi/2./prwt*(pbase-resp) 
c  CALCULATE  PROJECTILE  ACCELERATION 
254  z(l)=areab*(pbase-resp)/prwt 
if(z(l)  lt.0.0)go  to  257 
go  to  258 

257  if(iswl.eq.0)z(l)=0.0 

258  if(y(l).lL0.0)y(l)=0.0 
z(2)=y(l) 
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c  GET  BURNING  RATE 
do  264  m=l,nprop 
z(ibrp+m)=0.0 
d2xdt2(m)=0.0 
if(ibo(m).eq.l)  goto  264 
do  262  k=  1  ,nbr(m) 
if(pmean.gt.pres(m,k))go  to  262 
go  to  263 

262  continue 
k=nbr(m) 

263  pmix=pmean 

if(igrad.eq.  3)pmix=pbrch-(ak  1 1  *pbase+akl2)/6.*zb*zb 

if(igrad.eq.4)pmix=pbrch+(alt+a2t*pbase)*j3zb/vzb+bt*j4zbMb 

if(pmix.lL.99*pmean)pmix=pmean 

z(ibrp+m)=beta(m,k)*(pmix*  1  .e-6)**alpha(mjc) 

abr(m)=alpha(mjc) 

bbr(m)=beta(m,k) 

d2xdt2(m)=beta(mjc)*alpha(m,k)*(pmix*  1  .e*6) 

&  **(alpha(m4c)-l.)*dpmdt*l.e-6 

264  continue 

do  21  i=l,nde 
d(i)=(z(i)-b(j)*p(i))*a(j) 
y(i)=deltat*d(i)+y(i) 
p(i)=3.*d(i)-ak(j)*z(i)+p(i) 

21  continue 

1 1  continue 
t=t+deltat 
told=y(3) 

if(pmaxm.gtpmean)go  to  281 

pmaxm=pmean 

tpmaxm=y(3) 

281  if(pmaxba.gt.pbase)go  to  282 
pmaxba=pbase 
tpmxba=y(3) 

282  if(pmaxbr.gt.pbrch)go  to  283 
pmaxbr=pbrch 
tpmxbn=y(3) 

283  continue 
if(y(3).lt.ptime)go  to  272 
ptime=ptime-fdeltap 
pjt=y(2)+y(7) 

write(6,7)y(3),z(l),y(l),pjt,pmean,pbase,pbrch 

7  format(lx,7ell.4) 
if(igrad.gt.2)then 
pjt=y(2)+y(7) 
pit=y(10)+y(7) 
write(6,427)prt,pjt 

427  fonnatClx.’prop  traver,ell.4,’proj  travel’, el  1.4) 
endif 

272  continue 

if(t.gt.tstop)goto  200 
if(y(2)+y(7).gt.travp)go  to  200 
imvelo=y(l) 
tmvelo=y(3) 
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disto=y(2)+y(7) 
go  to  19 

200  write(6,311)t,y(3) 

311  format(lx,’  deltat  t’,  el4.6,  *  intg  t’,el4.6) 
write(6,3 1 2)pmaxm,tpmaxm 

312  format(lx,’PMAXMEAN  Pa  \el4.6,’  time  at  PMAXMEAN  sec  \el4.6) 
write(6,3 1 3)pmaxba,tpmxba 

313  format(  1  x, ’PM AXB ASE  Pa  \el4.6,’  time  at  PMAXBASE  sec  \el4.6) 
write(6,3 14)pmaxbr,tpmxbr 

314  format(  1  x.’PMAXBREECH  Pa  \el4.6,’  time  at  PMAXBREECH  sec  \el4.6) 
if(y(2)+y(7).le.travp)go  to  303 
dfract=(travp-disto)/(y(2)+y(7)-disto) 

rmvel=(y(  1  )-rmvelo)*dfract+rmvelo 
tmvel=(y(3)-tmvelo)*dfract+tmvelo 
write(6,3 1 8)rmvel,tmvel 

318  format(lx, ’muzzle  VELOCITY  m/s  ’,el4.6,’  time  of  muzzle  velocity  s 
&ec  ’,el4.6) 

goto  319 

303  write(6,327)y(l),y(3) 

327  format(lx, ’velocity  of  projectile  m/s  ’,el4.6,’  at  this  time  msec 
&’,el4.6) 

319  efi=chwi*forcig/(gamai-l.) 
efp=0.0 

do  315  i=  1  ,nprop 

efp=efp+chwp(i)’,'forcp(i)/(gamap(i)-l  .0) 

315  continue 
tenerg=efi+efp 
write(6,317)tenerg 

317  format(lx, ’total  initial  energy  available  J  =  ’,el4.6) 
tengas=chwi*forcig*tgas/(gamai- 1  .)Aempi 
do  135  i=l,nprop 

tengas=(frac(i)*chwp(i)*forcp(i)*tgas/tempp(i)/(gamap(i)- 1  .))+teng 
&as 

write(6,328)i,frac(i),tbo(i) 

328  foimat(’  FOR  PROPELLANT  ’,12,'  MASSFRACT  BURNT  IS  ’,el4.6 
&,’  at  time  in  sec  \el4.6) 

135  continue 
write(6,136)tengas 

136  foimat(lx, ’total  energy  remaining  in  gas  J=  ’,el4.6) 
write(6,320)elpt 

320  format(lx, ’energy  loss  from  projectile  translation  J=  ’,el4.6) 
write(6,321)elpr 

321  format(lx, ’energy  loss  from  projectile  rotation  J=  ’,el4.6) 
write(6,322)elgpm 

322  format(lx,’energy  lost  to  gas  and  propellant  motion  J=  ’,el4.6) 
write(6,323)elbr 

323  format( lx, ’energy  lost  to  bore  resistance  J=  ’,el4.6) 
write(6,324)elrc 

324  format(lx, ’energy  lost  to  recoil  J=  ’,el4.6) 
write(6,325)elht 

325  format(lx, ’energy  loss  from  heat  transfer  J=  ’,el4.6) 
write(6,326)elar 

326  format(lx,’energy  lost  to  air  resistance  J=  ’,el4.6) 
stop 
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20  write(*,140) 

140  format(lx,’end  of  file  encounter’) 
stop 

30  write(*, 150) 

999  continue 
998  continue 

150  format(lx,’read  or  write  error’) 
stop 
end 

subroutine  prf7 1 0(pd,gd,gl,np,x,frac,vol,surf1dsdx) 
common  nsljcpr,fracsl(10),dsdxsl(10),surfsl(10), 

&  nslp(10),tsl(10),pbrch,pbase,pmean,bbr(10),abr(10), 

&  deltat,yar(20)  'grad 
dimension  ts(10),coef(10) 
pi=3.141593 
nsl=0 
C 

C  pd=perforation  diameter 

C  gd=OUTER  DIA 

C  gl=GRAIN  LENGTH 

C  NP=NUMBER  OF  PERFS 

C 

C  SURF=OUTPUT  SURFACE  AREA 

C  frac=OUTPUT  MASS  FRACTION  OF  PROPELLANT  BURNED 
C 

C  w  =  web  =  distance  between  perforation  edges 
C  =  distance  between  outside  perf  edge  and  edge  of  grain 
C 

C  p  =  distance  between  perforation  centers 

C 

C  xl  =  distance  to  inner  sliver  burnout 

C 

C  x2  =  distance  to  outer  sliver  burnout  (frac=l) 

C 

if(np.eq.0)  go  to  2000 
IF(NP.EQ.l)GO  TO  3000 
IF(NF.eq.7)GO  TO  61 
if(np.eq.l9)go  to  4000 
if(np.eq.l5)go  to  5000 
60  WRITE(6,90) 

90  FORMAT(lX, ’UNACCEPTABLE  GRANULATION’) 

STOP 

61  w=(gd-3.*pd)/4. 

d=w+pd 
sqr3=sqrt(3.) 
xl=d/sqr3-pd/2. 
x2=(14.-3.*sqr3)*d/13.-pd/2. 
v0=pi/4  .*gl*(gd*gd-7*pd*pd) 
s0=2  *v0/gl+pi*gl*(gd+7  *pd) 
if  (x.gtw/2.+.0000001)  goto  20 
vol=pi/4.*(gl-2.*x)*((gd-2.*x)**2-7.*(pd+2.*x)**2) 
surf=2.*vol/(gl-2.*x)+pi*(gl-2.*x)*((gd-2.*x)+ 

&  7.*(pd+2.*x)) 
frac=l.-vol/v0 


ill 


dsdx=-4*pi*(gd+7.*pd-3.*gl+18.*x) 
dsdxsl(kpr>=dsdx 
fracsl(kpr)=frac 
surfsl(kpr)=surf 
return 
20  nsl=l 

coef(kpr)=0. 

if(igrad.eq.  1  .or.igrad.eq.2)go  to  726 
if(nslp(kpr).eq.l)goto  26 
tsl(kpr)=yar(3) 

ts(kpr)=w/2  .*(- 1  .+(pbrch/pniean)**abr(kpr))/ 

&  (bbr(kpr)*(pbase*  1  .e-6)**abr(kpr)) 

26  continue 

coef(kpr)=(ts(kpr)+tsl(kpr)-(deltat+yar(3)))As(kpr) 
if(coef(kpr).gt.  1  .)coef(kpr)=  1 . 
if(coef(kpr).lt.0.)coef(kpr)=0. 

726  if(x.ge.x2)goto  30 

sl=0. 
s2=0. 
vl=0. 
v2=0. 
dsldx=0. 
ds2dx=0. 

y=sqit((pd+2.*x)**2-d*d) 

theta=atan(y/d) 

al=theta/4.*(pd+2.*x)**2-d/4.*y 
if(x.ge.xl)goto  25 
v 1=374.  *(gl-2.*x) 

vl=vl*(2.*sqi3*d*d-pi*(pd+2.*x)**2+24.*al) 

sl=2.*vl/(gl-2.*x) 

sl=sl+3.*(gl-2.*x)*(pi-6.*theta)*(pd+2.*x) 

25  yl=sqrt((gd-2*x)**2-(5*d-2*(pd+2*x))**2) 

chi=atan(yl/(5.*d-2.*(pd+2.*x))) 
y2=sqrt((pd+2.*x)**2-(3.*d-2.*(pd+2.*x))**2) 
phi=atan(y2/(3*d-2.*(pd+2.*x») 
a2=phi*(pd+2.*x)**2-chi*(gd-2.*x)**2 
a2=(a2+2.*sqi3*d*sqrt((3*d-pd-2.*x)*(3*d-gd+2.*x)))/8. 
v2=pi*(gd-2.*x)**2-6.*sqi3*d*d4.*pi*(pd+2*x)**2 
v2=(v2+24.  *(a  1  +2.  *a2))*(gl-2.  *x)/4. 
s2=2.*v2/(gl-2.*x) 

s2=s2+(gl-2.*x)*((pi-6.*chi)*(gd-2*x)+2.*(2.*pi-3.*phi-3.*theta 
&  )*(pd+2.*x)) 
vol=vl+v2 
surf=sl+s2 
frac=l.-vol/v0 
dsdx=-surf/(x2-x) 

dsdx=coef(kpr)*dsdxsl(kpr)+(  1  ,-coef(kpr))*dsdx 
dsdxsl(kpr)=dsdx 

frac=coef(kpr)*fracsl(kpr)+(l  .-coef(kpr))",frac 
fracsl(kpr)=frac 

surf=coef(kpr)*surfsl(kpr)+(l  ,-coef(kpr))*surf 
surfsl(kpr)=surf 
return 
30  vol=0. 
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surf=0. 

frac=fracsl(kpr)*coef(kpr)+ 1  ,-coef(kpr) 

fracsl(kpr)=frac 

if(frac.gt..9999)  frac=l. 

if(frac.gt..9999)retum 

dsdx=0. 

dsdx=dsdxsl(kpr)*coef(kpr) 

dsdxsl(kpf)=dsdx 

if(abs(dsdx).lt.  1  .)dsdx=0. 

surf=sunsl(kpr)*coef(kpr) 

surfsl(kpr)=surf 

return 

C 

C  ZERO  PERF  CALCULATIONS  START  HERE. 

C 

2000  if(gd-2.*x.le.0.0)  go  to  2001 
v0=pi*gd*gd/4.*gl 
vol=pi*(gd-2.*x)**2/4.*(gl-2.*x) 
frac=l.-voVvO 

surf=pi/2*Cgd-2.*x)**2+pi*(gd-2*x)*(gl-2*x) 

dsdx='2.*pi*(gd+gl-6.*x) 

return 

2001  surf=0. 
frac=1.0 
vol=0. 
dsdx=0. 
nsl=l 
return 

c 

c  one  perf  calculation  starts  here 

c 

3000  if(gd-pd-4.*x.le.0.0)  goto  3001 
vfc=pi/4.*(gd*gd-pd*pd)*gl 
vol=pi/4.*((gd-2*x)**2-(pd+2*x)**2)*(gl-2*x) 
frac=l.-vol/vO 

surf=pi/2*((gd-2*x)**2-(pd+2*x)**2) 

surf=surf+pi*(gd-2.*x)*(gl-2.*x) 

surf=surf+pi*(pd+2.*x)*(gl-2.*x) 

dsdx=-4.*pi*(gd+pd) 

return 

3001  surf=0. 
frac=1.0 
vol=0. 
dsdx=0. 
nsl=l 
return 

C 

C  Below  is  the  calculation  for  the  cylindrical  19  perf  grain. 

C 

C  INPUT 

C 

C  P  =  PERF  DIAMETER 

C  D  =  GRAIN  DIAMETER 

C  GL  =  GRAIN  LENGTH 
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C  X  =  DISTANCE  BURNT 
C 

C  OUTPUT 
C 

C  VOL  =  THE  VOLUME  OF  ONE  GRAIN  AT  X. 

C  SURF  =  THE  SURFACE  AREA  OF  ONE  GRAIN  AT  X. 

C  FRAC  =  THE  FRACTION  OF  GRAIN  BURNT  AT  X. 

C 

C  W=WEB 
C 

4000  p=pd 
d=gd 

W=(D-5.*P)/6. 

PI=3. 14 1592654 
SQRT3=SQRT(3.) 

SQRT5=SQRT(5.) 

SQRT6=SQRT(6.) 

SQRT10=SQRT(10.) 

C 

C  INITIAL  VOLUME  AND  SURFACE  AREA 
C 

V0=PI/4.*GL*(D*D-19.*P*P) 

S0=2.*V0/GL+PI*GL*(D+19.*P) 

C 

C  XI  =  DISTANCE  TO  INNER  SLIVERR  BURNOUT 

C  X2  =  DISTANCE  TO  OUTER  SLIVER  BURNOUT 

C  DBC  =  DISTANCE  BETWEEN  PERFORATION  CENTERS 
C  ASSUMES  BURNOUT  DOES  NOT  OCCUR  IN  LONGITUDINAL  DIRECTION 
C  W1  =  SECONDARY  WEB 
C 

DBC=W+P 

W1=0.5*(D-P-2.*SQRT3*DBC) 

X 1 =DBC/SQRT3-P/2. 

X2=0.25*(DBC*(6.-SQRT10)-2.*P) 

IF(X.GT.W/2.)GO  TO  110 
C 

C  NOT  SLIVERED  YET 
C 

VOL=PI/4.*(GL-2.*X)*((D-2*X)**2-19*(P+2*X)**2) 

SURF=2.*VOL/(GL-2.*X)+PI*(GL-2.*X)*(D-2.*X+19.*(P+2.*X)) 

dsdx=pi*(-4*D+36*GL-76*P-216*x) 

FRAC=  1  .-VOL/VO 
dsdxsl(kpr)=dsdx 
firacsl(kpr)=frac 
surfslOcpr)=surf 
RETURN 
C 

C  Vl=TOTAL  VOLUME  OF  INNER  SLIVER,  V2=TOTAL  VOLUME  OF  OUTER  SLIVER 
C  Sl=TOTAL  SURFACE  AREA  OF  INNER  SLIVERS,  S2=TOTAL  SURFACE  AREA  OF 
C  OUTER  SLIVERS 


S2=0. 

DELTA=0. 

CHI=0. 

NSL=1 

coef(kpr)=0. 

if(igrad.eq.l.or.igrad.eq.2)go  to  111 
if(nslp(kpr).eq.l)goto  728 
tsl(lpr)=yar(3) 

ts(kpr)=w/2  .*(- 1  ,+(pbrch/pmean)**abr(kpr))/ 

&  (bbr(kpr)*(pbase*  1  .e-6)**abr(kpr)) 

728  continue 

coef(kpr)=(ts(kpr)+tsl(kpr)-(deltat+yar(3)))As(kpr) 
if(coef(kpr).gt.  1  .)coef(kpr)=  1 . 
if(coef(kpr).lt.0.)cocf(kpr>=0. 

727  A3=0. 

IF(X.GE.X2)GO  TO  130 
THETA=ACOS(DBC/(P+2.*X)) 

A1=THETA/4.*(P+2.*X)**2-DBC/4.*SQRT((P+2.*X)**2-DBC*DBC) 
IF(X.GT.Xl)GO  TO  120 

V1=3*(GL-2*X)*(2*SQRT3*DBC*DBC-PI*(P+2.*X)**2+24*A1) 

S1=2*V1/(GL-2*X)+12.*(GL-2.*X)*(PI-6.*THETA)*(P+2*X) 

120  PHI=ACOS((5.*D-13*P-36*X)/(12*(P+2.*X))) 

XI=ACOS((13*D-5.*P-36.*X)/(12.*(D-2.*X))) 

IF(X.LE.Wl/2.)GO  TO  125 
DELTA=ACOS((2.*D-P-6.*X)/(SQRT3*(D-2.*X))) 
CHI=ACOS((D-2.*P-6.*X)/(SQRT3*(P+2*X))) 
A3=.125*(CHI*(P+2.*X)**2-DELTA*(D-2.*X)**2 
&  +2.*SQRT6*DBC*SQRT(6.*DBC*(P+2*X-DBC)-(P+2*X)**2)) 
125  A2=.125*(PHI*(P+2*X)**2-XI*(D-2.*X)**2 

&  +2.*SQRT5*DBC*SQRT((5*DBC-P-2.*X)*(5*DBC-D+2*X))) 
V2=.25*(GL-2.*X)*(PI*(D-2.*X)**2-7*PI*(P+2.*X)**2 
&  -24.*SQRT3*DBC*DBC+48.*(A1+A2+A3)) 
S2=2.*V2/(GL-2.*X)-KGL-2.*X)*((D-2.*X)*(PI-6.*(XI+DELTA» 

&  +(P+2.*X)*(7.*PI-6.*(2.*THETA+CHI+PHI))) 

VOL=Vl+V2 

SURF=S1+S2 

DSDX=-SURF/(X2-X) 

FRAC=1. -VOL/VO 

dsdx=coef(kpr)*dsdxsl(kpr)+(  1  .-coef(kpr))*dsdx 
dsdxsl(kpr)=dsdx 

firac=coef(kpr)*fracsl(kpr)+(  1  .-coef(kpr))*frac 
fracsl(kpr)=frac 

surf=coef(kpr)*surfsl(kpr)+(  1  .-coef(kpr))*surf 
surfsl(kpr)=surf 
RETURN 
130  VOL=0. 

SURF=0. 

frac=fracsl(kpr)*coef(kpr)+ 1  ,-coef(kpr) 

fracsl(kpr)=frac 

if(frac.gt..9999)  frac=l. 

if(frac.gt..9999)retum 

dsdx=0. 

dsdx=dsdxsl(kpr)*coef(kpr) 

dsdxsl(kpr)=dsdx 
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on  non  no  n  no  on  n  non  nnnnn 


if(abs(dsdx).lt  1  .)dsdx=0. 
surf=surfsl(kpr)*coef(kpr) 
surfsl(kpr)=surf 
RETURN 
C 

C  Below  is  the  calculation  for  the  19  perf  hex  grain. 

C 

C 

C  Translation  of  the  input  values. 

C  p=  perf  diameter 

C  d=  grain  diameter 

C  gl=  grain  length 

C  x=  distance  burnt 

C 

C  Translation  of  the  output  values. 

C  vol=  volume  of  one  grain  at  x. 

C  surf=  surface  area  of  one  grain  at  x. 

C  frac=  mass  fraction  of  the  grain  burnt  at  x. 

C 

C  Assignment  statement  for  pi. 

5000  pi=3. 141592654 
sqrt3=sqrt(3.) 

p=pd 
d=gd 


d=6w  +  5p  is  the  statement  for  die  grain  diameter  which  will  be 
used  to  calculate  the  web. 

To  calculate  the  web. 
w=  (d-5.*p)/6. 


Below  is  the  equation  to  calculate  the  distance  between  the  perf  cen¬ 
ters. 

dpc=  p  +  w 

To  calculate  the  grain  diameter  between  the  flats. 
f=  2.*(sqrt3*dpc  +  p/2.  +  w) 

To  calculate  the  distance  burnt 
Xl=dpc/sqrt3-p/2. 

X2=  0.125*(5.*dpc-4.*p) 

To  calculate  the  area. 

A=sqrt3/3.*((w+p/2.)**2)-pi/6.*((w+p/2.)**2) 

To  calculate  the  initial  volume  of  the  sharp  comer  grain. 
Vs=gl/4.*(2.*sqrt3*f**2-19.*pi*p**2) 

To  calculate  the  volume  that  will  be  removed  from  the  grain. 
Vi*6.*A*gl 

To  calculate  the  initial  volume  for  the  19hex  grain  with  rounded 
comers. 

Vo=  Vs  -  Vr 

To  calculate  the  initial  surface  area  of  the  sharp  comer  grain. 


no  o  non  nn  n  n  nn  non  nn 


Ss=2 .  *  V  s/gl+gl*(2 .  *sqrt3*f+ 1 9.  *pi*p) 


To  calculate  the  surface  area  that  will  be  removed  from  the  grain. 
Sr=  1 2.  *  A+gl*(w+p/2.)*(4  *sqrt3-2.  *pi) 

To  calculate  the  initial  surface  area  for  the  19hex  grian  with  rounded 
comers. 

So=  Ss-Sr 


To  calculate  the  unknows  of  the  grain  under  the  condition  x.le.5*w. 
if(0.1e.x.and.x.le.w/2.)  then 
A=sqrt3/3.*(w-2.*x+(p+2.*x)/2.)**2-pi/6.* 

&(w-2.*x+(p+2.*x)/2.)**2 

To  calculate  the  volume  that  will  be  removed  from  the  sharp  comer  grain. 
Vn=6.*A*(gl-2*x) 

To  calculate  the  volume  for  the  sharp  comer  grain  at  some  distance  burnt. 
Vn=.25*(gl-2.*x)*(2.*sqrt3*(f-2.*x)**2. 

&  -19.*pi*(pt-2.*x)**2.) 

To  calculate  the  volume  for  the  19hex  grain  with  rounded  comers. 

V=  Vn-Vr 


To  calculate  the  surface  area  that  will  be  removed  from  the  sharp 
comer  grain. 

Sn=12.*A+(gl-2.*x)*(w-2*x+(p+2.*x)/2.)*(4.*sqrt3-2.*pi) 

To  calculate  the  surface  area  for  the  shaip  comer  grain. 
Sn=2.*V/(gl-2.*x)+(gl-2.*x)*(2.*sqrt3*(f-2.*x)+ 

&  19.*pi*(p+2.*x)) 

To  calculate  the  surface  area  for  19hex  grain  with  rounded  comers. 

S=  Sn-Sr 
C  To  calculate  the  mass  fraction. 
frac=  1-V/Vo 

dsdx=-8.*sqrt3*(f-2.*x)-76.*pi*(p+2.*x)+(gl-2.*x)*(-4.*sqrt3+38. 

&  *pi)+16*sqrt3*(w+p/2.-x)-8*pi*(w+p/2.-x)+(gl-2.*x)* 

&  (4.*sqrt3-2.*pi) 
surf=S 
vol=V 

dsdxsl(kpr)=dsdx 

firacsl(kpr)=frac 

surfsl(kpr)=surf 

return 

endif 

C 

C  Due  to  the  cross  section  at  the  sliver  point  x=.5*w  there  will  be  24 
C  identical  inner  slivers,  12  identical  side  slivers.  After  slivering  the 
C  surface  area  and  the  volume  function  become  more  complex.  Each  type  of 
C  sliver  will  be  treated  seperately  and  later  the  volumes  will  be  combined 
C  to  complete  the  function. 

C 

C  To  calculate  the  12  identical  side  slivers  for  the  grain  x=.5/w. 
nsl=l 

coef(kpr)=0. 

if(igrad.eq.l.or.igrad.eq.2)go  to  729 
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if(nslp(kpr).eq.l)goto  730 
tsl(kpr)=yar(3) 

ts(kpr)=w/2  *(- 1  .+(pbrch/pmean)**abr(kpr))/ 

&  (bbr0cpr)*(pbase*l.e-6)**abr(kpr)) 

730  continue 

coef(kpr)=(ts(kpr)+tsl(kpr)-(deltat+yar(3)))/ts(kpr) 
if(coef(kpr).gt.  1  .)coef(kpr)=  1 . 
if(coef(kpr).lt.0.)coef(kpr)=0. 

729  if(w/2.1t.x.and.x.lt.Xl.and.x.lt.X2)  then 
C 

C  To  calculate  the  areas  of  the  grain. 

A=sqrt3/3.*(w-2.*x+(p+2.*x)/2.)**2-pi/6.* 

&(w-2.*x+(p+2.*x)/2.)**2 

theta=acos(dpc/(p+2.*x)) 

Al=theta/4.*(p+2.*x)**2-dpc/4.*sqrt((p+2.*x)**2-dpc**2) 
omega=acos(2.*dpc/(p+2.*x)- 1 .) 

A2=0.125*(p+2  *x)*((p+2  *x)*(omega+sin(omega))-2  *dpc*sin(omega» 
C  To  calculate  the  volumes  of  the  grain. 

Vl=3.*(gl-2.*x)*(2.*sqrt3*dpc**2-pi*(p+2.*x)**2+24.*Al) 

V2=6.*(gl-2.*x)*(2.*dpc**2-dpc*(p+2.*x)-pi/4.*(p+2.*x)**2 

&+2.*Al44.*A2) 

C  To  calculate  the  surface  areas  of  the  grain. 

Sl=2.*Vl/(gl-2.*x)+12.*(gl-2.*x)!it(pi-6.*theta)*(p+2.*x) 

S2=2.*V2/(gl-2.*x)+12.*(gl-2.*x)*(dpc+(p+2.*x)*(pi/2.-omega 

&-theta-sin(omega))) 

C  To  calculate  the  total  volume  and  total  surface  area. 

Vf=Vl+V2 

Sf=Sl+S2 

C  To  calculate  the  mass  fraction. 
fiac=1.-Vf/Vo 
surf=Sf 

dsdx=-surf/(x2-x) 

vol=Vf 

dsdx=coef(kpr)*dsdxsl(kpr)+(  1  .-coef(kpr))*dsdx 
dsdxsl(kpr)=dsdx 

firac=coef(kpr)*fracsl(kpr)+(l  .-coef(kpr))*frac 
fracsl(kpr)=frac 

surf=coef(kpr)*surfsl(kpr)+(l.-coef(kpr))*surf 

surfsl(kpr)=surf 

return 

endif 

if(x.gtXl  .and.x.lt.X2)then 
C  To  calculate  the  area  of  the  grain. 

A=sqrt3/3.*(w-2.*x+(p+2.*x)/2.)**2-pi/6.* 

&(w-2.*x+(p+2.*x)/2.)**2 

theta=acos(dpc/(p+2.*x)) 

Al=theta/4.*(p+2.*x)**2-dpc/4.*sqrt((p+2.*x)**2-dpc**2) 

omega=acos(2.*dpc/(p+2.*x)-l.) 

A2=0.125*(p+2.*x)*((p+2.*x)*(omega+sin(omega))-2.*dpc*sin(omega)) 
C  To  calculate  the  volume  of  the  grain. 

V2=6.*(gl-2.*x)*(2.*dpc**2-dpc*(p+2.*x)-pi/4.*(p+2.*x)**2 

&+2.*Al+4.*A2) 

C  To  calculate  the  surface  area  of  the  grain. 

S2=2.*V2/(gl-2.*x)+12.*(gl-2.*x)*(dpc+(p+2.*x)*(pi/2.-omega 


118 


&-theta-sin(omega))) 

C  To  calculate  the  volume  and  the  surface  area. 
Vf=V2 
Sf=S2 

C  To  calculate  the  the  mass  fraction. 
frac=l-Vf/Vo 
surf-Sf 

dsdx=-surf/(x2-x) 

vol=Vf 

dsdx=cocf(kpr)*dsdxsl(kpr)+(  1  .-coef(kpr))*dsdx 
dsdxsl(kpr)=dsdx 

frac=coef(kpr)*fracsl(kpr)-K  1  .-coef(kpr))*frac 
fracsl(kpr)=frac 

surf=coef(kpr)*surfsl(kpr)+(l.-coef(kpr))*surf 

surfsl(kpr)=surf 

return 

endif 

if(x.gtX2)then 

dsdx=0. 

surf=0. 

vol=0. 

firac=fracsl(kpr)*coef(kpr)+l  .-coef(kpr) 

firacsl(kpr)=frac 

if(frac.gL.9999)  frac=l. 

if(frac.gt..9999)retum 

dsdx=0. 

dsdx=dsdxsl(kpr)*coef(kpr) 
dsdxsl(kpr)=dsdx 
if(abs(dsdx).lt  1  .)dsdx=0. 
surf=surfsl(kpr)*coef(kpr) 
surfsl(kpr)=surf 
return 
endif 
end 
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Intentionally  left  blank. 
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APPENDIX  6: 

INPUT  DATA  FOR  IBRGA  TEST 
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Intentionally  left  blank. 
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9832.2384  12.7  12.7  1.0  0.0  457.2  4  1.00 
3  0.  15.381732  46.482  15.381732  54.102  12.7 
9.796  0  0.0  0.0 

5  0.0  0.0  0.0  .6  0.0  1.3  0.0  300.  0.  457. 

I.e20  2  3.0e+4  0.0  8.0e+5  0.2 
.001135  .01143  .46028  273.  1.  7.8612 
106.34  .9755  294.  .004712  1.4 

1 

1135.84  3142.75  .97544  9.7959  1.6605  1.23  7  3.175  .0508  1.23063 
1  1.0  .1105187  689.476 
.01  .1  30. 
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Intentionally  left  blank. 
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APPENDIX  7: 
OUTPUT  FROM  IBRGA 
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Intentionally  left  blank. 
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the  input  file  is  i6315a 
using  rga  gradient 


chamber  distance  cm  chamber  diameter  cm 


0.153817E+02 

0.153817E+02 

0.127000E+02 

0.982089E+04 


O.OOOOOOE+OO 
0.464820E+02 
0.541020E+02 
chamber  volume  cm**3 
groove  diam  cm  0.127000E+02 
land  diam  cm  0.127000E+02 
groove/land  ratio  0.100000E+01 
twist  tums/caliber  0.000000E+00 


projectile  travel  cm  0.457200E+03 
gradient  #  4 

friction  factor  0.100000E+01 


projectile  mass  kg  0.979600E+01 
switch  to  calculate  energy  lost  to  air  resistance  J  0 
fraction  of  woik  against  bore  used  to  heat  the  tube  O.OOOOOOE+OO 
gas  pressure  Pa  O.OOOOOOE+OO 
number  barrel  resistance  points  5 
bore  resistance  MPa  -  travel  cm 
O.OOOOOOE+OO  O.OOOOOOE+OO 
O.OOOOOOE+OO  0.600000E+00 
O.OOOOOOE+OO  0.130000E+01 
O.OOOOOOE+OO  0.300000E+03 
O.OOOOOOE+OO  0.457000E+03 

mass  of  recoiling  parts  kg  0.100000E+21 
numoer  of  recoil  point  pairs  2 
recoil  force  N  recoil  time  sec 
0.300000E+05  O.OOOOOOE+OO 

0.800000E+06  0.200000E+00 

free  convective  heat  transfer  coefficient  w/cm**2  K  0.113500E-02 
chamber  wall  thickness  cm  0.114300E-01 
heat  capacity  of  steel  of  chamber  wall  J/g  K  0.460280E+00 
initial  temperature  of  chamber  wall  K  0.273000E+03 
heat  loss  coefficient  0.100000E+01 
density  of  chamber  wall  steel  g/cm**3  0.786120E+01 


impetus  of  igniter  propellant  J/g  0.106340E+03 
covolume  of  igniter  cm**3/g  0.975500E+00 
adiabatic  flame  temperature  of  igniter  propellant  K  0.294000E+03 
initial  mass  of  igniter  kg  0.471200E-02 
ratio  of  specific  heats  for  igniter  0.140000E+01 
there  are  1  propellants 
for  propellant  number  1 
impetus  of  propellant  J/g  0.113584E+O4 
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adiabatic  temperature  of  propellant  K  0.314275E+04 

covolume  of  propellant  cm**3/g  0.975440E+00 

initial  mass  of  propellant  kg  0.979590E+01 

density  of  propellant  g/cm**3  0.166050E+01 

ratio  of  specific  heats  for  propellant  0.123000E+01 

number  of  perforations  of  propellant  7 

length  of  propellant  grain  cm  0.317500E+01 

diameter  of  perforation  in  propellant  grains  cm  0.508000E-01 

outside  diameter  of  propellant  grain  cm  0.123063E+01 


for  propellant  1  the  total  number  of  grains  is  0.158099E+04 
number  of  burning  rate  points  1 

exponent  coefficient  pressure 

cm/scc-MPa**ai  MPa 

0.100000E+01  0.110519E+00  0.689476E+03 

time  increment  msec  0.100000E-01  print  increment  msec  0.100000E+00 

time  to  stop  calculation  msec  0.300000E+02 

area  bore  mA2  0.126677E*01 

pressure  from  ign  pa  0.127925E+06 

volume  of  unbumt  prop  mA3  0.589937E-02 

init  cham  vol-cov  ign  mA3  0.981629E-02 

time  acc  vel  dis  mpress  pbase  pbrch 
0.1000E-04  0.1696E+03  0.1399E-02  0.5582E-08  0.1312E+06  0.1312E+06  0.1313E+06 
prop  travel  0.8638E-llproj  travel  0.5582E-08 

0.1100E-03  0.2182E+03  0.2070E-01  0.1070E-05  0.1688E+06  0.1688E+06  0.1689E+06 
prop  travel  0.1936E-08proj  travel  0.1070E-05 

0.2100E-03  0.2790E+03  0.4545E-01  0.4327E-05  0.2159E+06  0.2158E+06  0.2160E+06 
prop  travel  0.8074E-08proj  travel  0.4327E-05 

0.3100E-03  0.3543E+03  0.7699E-01  0.1039E-04  0.2742E+06  0.2740E+06  0.2743E+06 
prop  travel  0.2001 E-07proj  travel  0.1039E-04 

0.4100E-03  0.4464E+03  0.1169E+00  0.2000E-04  0.3454E+06  0.3452E+06  0.3455E+06 
prop  travel  0.4001 E-07proj  travel  0.2000E-04 

0.5100E-03  0.5576E+03  0.1669E+00  0.3410E-04  0.4315E+06  0.4312E+06  0.4317E+06 
prop  travel  0.7140E-07proj  travel  0.3410E-04 

0.6100E-03  0.6904E+03  0.2291E+00  0.5379E-04  0.5344E+06  0.5339E+06  0.5346E+06 
prop  travel  0.1190E-06proj  travel  0.5379E-04 

0.7100E-03  0.8474E+03  0.3058E+00  0.8040E-04  0.6560E+06  0.6553E+06  0.6563E+06 
prop  travel  0.1898E-06proj  travel  0.8040E-04 

0.8100E-03  0.1031E+04  0.3995E+00  0.1155E-03  0.7986E+06  0.7975E+06  0.7990E+06 
prop  travel  0.2943E-06proj  travel  0.1155E-03 

0.9100E-03  0.1245E+04  0.5130E+00  0.1610E-03  0.9646E+06  0.9631E+06  0.9653E+06 
prop  travel  0.4479E-06proj  travel  0.1610E-03 

0.1010E-02  0.1493E+04  0.6497E+00  0.2189E-03  0.1157E+07  0.1155E+07  0.1158E+07 
prop  travel  0.6734E-06proj  travel  0.2189E-03 

0.1110E-02  0.1780E+04  0.8130E+00  0.2918E-03  0.1380E+07  0.1376E+07  0.1381E+07 
prop  travel  0.1005E-05proj  travel  0.291 8E-03 

0.1210E-02  0.2110E+04  0.1007E+01  0.3825E-03  0.1637E+07  0.1632E+07  0.1639E+07 
prop  travel  0.1494E-05proj  travel  0.3825E-03 

0.1310E-02  0.2492E+04  0.1237E+01  0.4944E-03  0.1934E*-07  0.1927E+07  0.1936E+07 
prop  travel  0.2217E-05proj  travel  0.4944E-03 

0.1410E-02  0.293 1E+04  0.1507E+01  0.6312E-03  0.2276E+07  0.2266E+07  0.2281E+07 
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prop  travel  0.3292E-05proj  travel  0.6312E-03 

0.1510E-02  0.3437E+04  0.1825E+01  0.7974E-03  0.2673E+07  0.2658E+07  0.2679E+07 
prop  travel  0.4893E-05proj  travel  0.7974E-03 

0.1610E-02  0.402 1E+04  0.2197E+01  0.9981E-03  0.3131E+07  0.3110E+07  0.3141E+07 
prop  travel  0.7289E-05proj  travel  0.998 IE-03 

0.1710E-02  0.4694E+04  0.2632E+01  0.1239E-02  0.3662E+07  0.3630E+07  0.3676E+07 
prop  travel  0.1089E-04proj  travel  0.1239E-02 

0.1810E-02  0.5470E+04  0.3140E+01  0.1527E-02  0.4277E+07  0.4230E+07  0.4298E+07 
prop  travel  0.1630E-04proj  travel  0.1527E-02 

0.1910E-02  0.636 1E+04  0.3730E+01  0.1870E-02  0.4989E+07  0.4919E+07  0.5020E+07 
prop  travel  0.2447E-04proj  travel  0.1870E-02 

0.2010E-02  0.7384E+04  0.4416E+01  0.2276E-02  0.5814E+07  0.5710E+07  0.5859E+07 
prop  travel  0.3681E-04proj  travel  0.2276E-02 

0.2110E-02  0.8556E+04  0.5212E+01  0.2757E-02  0.6769E+07  0.6616E+07  0.6836E+07 
prop  travel  0.5544E-04proj  travel  0.2757E-02 

0.2210E-02  0.9892E+04  0.6133E+01  0.3323E-02  0.7874E+07  0.7650E+07  0.7974E+07 
prop  travel  0.8353E-04proj  travel  0.3323E-02 

0.2310E-02  0.U41E+05  0.7196E+01  0.3988E  02  0.9153E+07  0.8825E+07  0.9300E+07 
prop  travel  0.1258F-03proj  travel  0.3988E-02 

0.2410E-02  0.1313E+05  0.8421E+01  0.4767E-02  0.1063E+08  0.1015E+08  0.1085E+08 
prop  travel  0.1891E-03proj  travel  0.4767E-02 

0.2510E-02  0.1507E+05  0.9829E+01  0.5678E-02  0.1234E+08  0.1165E+08  0.1265E+08 
prop  travel  0.2835E-03proj  travel  0.5678E-02 

0.2610E-02  0.1724E+05  0.1144E+02  0.6740E-02  0.1431E+08  0.1333E+08  0.1475E+08 
prop  travel  0.4232E-03proj  travel  0.6740E-02 

0.2710E-02  0. 1967E+05  0.1328E+02  0.7974E-02  0.1658E+08  0.1521E+08  0.1719E+08 
prop  travel  0.6278E-03proj  travel  0.7974E-02 

0.2810E-02  0.2238E+05  0.1538E+02  0.9405E-02  0.1918E+08  0.1731E-K)8  0.2003E+08 
prop  travel  0.9242E-03proj  travel  0.9405E-02 

0.2910E-02  0.254 1E+05  0.1777E+02  0.1106E-01  0.2217E+08  0.1965E+08  0.2332E+08 
prop  travel  0.1348E-02proj  travel  0.1106E-01 

0.3010E-02  0.2882E+05  0.2047E+02  0.1297E-01  0.2560E+08  0.2229E+08  0.2710E+08 
prop  travel  0.1943E-02proj  travel  0.1297E-01 

0.3110E-02  0.3268E+05  0.2354E+02  0.1517E-01  0.2951E+08  0.2527E+08  0.3144E+08 
prop  travel  0.2765E-02proj  travel  0.1517E-01 

0.3210E-02  0.3709E+05  0.2702E+02  0.1769E-01  0.3397E+08  0.2869E+08  0.3637E+08 
prop  travel  0.3877E-02proj  travel  0.1769E-01 

0.3310E-02  0.4217E+05  0.3098E+02  0.2059E-01  0.3903E+08  0.3261E+08  0.4194E+08 
prop  travel  0.535 1E-O2proj  travel  0.2059E-01 

0.3410E-02  0.4802E+05  0.3548E+02  0.2391E-01  0.4475E+08  0.3714E+08  0.4820E+08 
prop  travel  0.7263E-02proj  travel  0.2391E-01 

0.3510E-02  0.5475E+05  0.4060E+02  0.2770E-01  0.5119E+08  0.4234E+08  0.5518E+08 
prop  travel  0.9691E-02proj  travel  0.2770E-01 

0.3610E-02  0.6243E+05  0.4644E+02  0.3205E-01  0.5841E+08  0.4828E+08  0.6294E+08 
prop  travel  0.1271E-01proj  travel  0.3205E-01 

0.3710E-02  0.7109E+05  0.5310E+02  P.3702E-01  0.6645E+08  0.5497E+08  0.7153E+08 
prop  travel  0.1641E-01proj  travel  0.3702E-01 

0.3810E-02  0.8072E+05  0.6068E+02  0.4270E-01  0.7536E+08  0.6242E+08  0.8100E+08 
prop  travel  0.2085E-01proj  travel  0.4270E-01 

0.3910E-02  0.9131E+05  0.6926E+02  0.4919E-01  0.8514E+08  0.7061E+08  0.9139E+08 
prop  travel  0.261  lE-Olproj  travel  0.4919E-01 

0.401 0E-02  0.1028E+06  0.7895E-1-02  0.5659E-01  0.9580E+08  0.7948E+08  0.1027E+09 
prop  travel  0.3228E-01proj  travel  0.5659E-01 

0.4110E-02  0.1151E+06  0.8983E+02  0.6502E-01  0.1073E+09  0.8898E+08  0.1150E+09 


129 


prop  travel  0.3944E-01proj  travel  0.6502E-01 

0.4210E-02  0.1280E+06  0.1020E+03  0.7460E-01  0.1197E+09  0.9901E408  0.1282E+09 
prop  travel  0.4767E-01proj  travel  0.7460E-01 

0.4310E-02  0.1416E+06  0.1154E+03  0.8546E-01  0.1327E+09  0.1095E+09  0.1421E+09 
prop  travel  0.5707E-01proj  travel  0.8546E-01 

0.4410E-02  0.1555E+06  0.1303E+03  0.9773E-01  0.1464E+09  0.1203E+09  0.1568E+09 
prop  travel  0.6774E-01proj  travel  0.9773E-01 

0.4510E-02  0.1697E+06  0.1465E+03  0.1116E+00  0.1605E+09  0.1313E+09  0.1721E+09 
prop  travel  0.7980E-01proj  travel  0.1116E+00 

0.4610E-02  0.1840E+06  0.1642E+03  0.1271E+00  0.1748E+09  0.1423E+09  0.1877E+09 
prop  travel  0.9336E-01proj  travel  0.1271E+00 

0.4710E-02  0.1981E+06  0.1833E+03  0.1444E+00  0.1893E+09  0.1532E+09  0.2035E+09 
prop  travel  0.1085E+00proj  travel  0.1444E+00 

0.4810E-02  0.2118E+06  0.2038E+03  0.1638E+00  0.2036E+09  0.1638E+09  0.2192E+09 
prop  travel  0.1254E+00proj  travel  0.1638E+00 

0.4910E-02  0.2249E+06  0.2256E+03  0.1852E+00  0.2175E+09  0.1739E+09  0.2347E+09 
prop  travel  0.1442E+00proj  travel  0.1852E+00 

0.5010E-02  0.2372E+06  0.2487E+03  0.2089E+00  0.2308E+09  0.1834E+09  0.2496E+09 
prop  travel  0.1648E+00proj  travel  0.2089E+00 

0.5110E-02  0.2485E+06  0.2730E+03  0.2350E+00  0.2433E+09  0.1922E+09  0.2637E+09 
prop  travel  0.1876E+00proj  travel  0.2350E+00 

0.5210E-02  0.2587E+06  0.2983E+03  0.2636E+00  0.2547E+09  0.2001E+09  0.2769E+09 
prop  travel  0.2124E+00proj  travel  0.2636E+00 

0.5310E-02  0.2677E+06  0.3247E+03  0.2947E+00  0.2651E+09  0.2070E+09  0.2889E+09 
prop  travel  0.2395E+00proj  travel  0.2947E+00 

0.5410E-02  0.2754E+06  0.3518E+03  0.3285E+00  0.2741E+09  0.2129E+09  0.2998E+09 
prop  travel  0.2689E+00proj  travel  0.3285E+00 

0.5510E-02  0.2817E+06  0.3797E+03  0.3651E+00  0.2819E+09  0.2178E+09  0.3092E+09 
prop  travel  0.3006E+00proj  travel  0.3651E+00 

0.5610E-02  0.2867E+06  0.4081E+03  0.4045E+00  0.2883E+09  0.2217E+09  0.3172E+09 
prop  travel  0.3348E+00proj  travel  0.4045E+00 

0.5710E-02  0.2903E+06  0.4370E+03  0.4467E+00  0.2933E+09  0.2245E+09  0.3239E+09 
prop  travel  0.3714E+00proj  travel  0.4467E+00 

0.5810E-02  0.2928E+06  0.4661E+03  0.4919E+00  0.2970E+09  0.2264E+09  0.3291E+09 
prop  travel  0.4105E+00proj  travel  0.4919E+00 

0.5910E-02  0.2940E+06  0.4955E+03  0.5400E+00  0.2995E+09  0.2274E+09  0.3329E+09 
prop  travel  0.4520E+00proj  travel  0.5400E+00 

0.6010E-02  0.2942E+06  0.5249E+03  0.5910E+00  0.3008E+09  0.2275E+09  0.3355E+09 
prop  travel  0.496 lE+OOproj  travel  0.5910E+00 

0.6110E-02  0.2934E+06  0.5543E+03  0.6449E+00  0.3010E+09  0.2269E+09  0.3369E+09 
prop  travel  0.5428E+00proj  travel  0.6449E+00 

0.6210E-02  0.2918E+06  0.5835E+03  0.7018E400  0.3003E+09  0.2257E+09  0.3371E+09 
prop  travel  0.5919E+00proj  travel  0.7018E+00 

0.6310E-02  0.2894E+06  0.6126E+03  0.7616E+00  0.2988E+09  0.2238E+09  0.3364E+09 
prop  travel  0.6435E+00proj  travel  0.7616E+00 

0.6410E-02  0.2864E+06  0.6414E+03  0.8243E+00  0.2965E+09  0.2215E+09  0.3348E+09 
prop  travel  0.6976E+00proj  travel  0.8243E+00 

0.6510E-02  0.2828E+06  0.6699E+03  0.8899E+00  0.2936E+09  0.2187E+09  0.3323E+09 
prop  travel  0.7542E+00proj  travel  0.8899E+00 

0.6610E-02  0.2788E+06  0.6979E+03  0.9583E+00  0.2901E+09  0.2156E+09  0.3293E+09 
prop  travel  0.8132E+00proj  travel  0.9583E+00 

0.6710E-02  0.2744E+06  0.7256E+03  0.1029E401  0.2862E+09  0.2122E+09  0.3256E+09 
prop  travel  0.8746E+00proj  travel  0.1029E+01 

0.6810E-02  0.2697E+06  0.7528E+03  0.1103E401  0.2820E+09  0.2085E+09  0.3214E+09 


130 


prop  travel  0.9384E+0Oproj  travel  0.1103E+01 

0.6910E-02  0.2647E+06  0.7795E+03  O.II8OE4OI  0.2774E+09  0.2047E+09  0.3169E+09 
prop  travel  0.1005E+01proj  travel  0.1180E+01 

0.7010E-02  0.2596E+06  0.8058E+03  0.1259E+O1  0.2727E+09  0.2008E+09  0.3120E+09 
prop  travel  0.1073E+01proj  travel  0.1259E+01 

0.7110E-02  0.2544E+06  0.8315E+03  0.1341E+01  0.2677E+09  0.1968E+09  0.3069E+09 
prop  travel  0.1144E+01proj  travel  0.1341E+01 

0.7210E-02  0.2492E+06  0.8566E+03  0.1426E+01  0.2627E+09  0.1927E+09  0.3016E+09 
prop  travel  0.1217E+01proj  travel  0.1426E+01 

0.7310E-02  0.2439E+06  0.8813E+03  0.1512E+01  0.2576E+09  0.1886E+09  0.2961E+09 
prop  travel  0.1292E+01proj  travel  0.1512E+01 

0.7410E-02  0.2385E+06  0.9054E+03  0.1602E+01  0.2525E+09  0.1845E+09  0.2905E+09 
prop  travel  0.1369E+01proj  travel  0.16Q2E+01 

0.7510E-02  0.2333E+06  0.9290E+03  0.1694E+01  0.2473E+09  0.1804E+09  0.2849E+09 
prop  travel  0.1449E+01proj  travel  0.1694E+01 

0.7610E-02  0.228 1E+06  0.9521E+03  0.1788E+01  0.2422E+09  0.1764E+09  0.2793E+09 
prop  travel  0.1530E+01proj  travel  0.1788E+01 

0.7710E-02  0.2229E+06  0.9746E+03  0.1884E+01  0.2372E+09  0.1724E+09  0.2737E+09 
prop  travel  0.1614E+01proj  travel  0.1884E+01 

0.7810E-02  0.2179E+06  0.9967E+03  0.1983E+01  0.2322E+09  0.1685E+09  0.2681E+09 
prop  travel  0.1700E+01proj  travel  0.1983E+01 

0.7910E-02  0.2129E+06  0.1018E+04  0.2083E+01  0.2272E+09  0.1646E+09  0.2625E+09 
prop  travel  0.1787E+01proj  travel  0.2083E+01 

0.8010E-02  0.208 1E+06  0.1039E+04  0.2186E+01  0.2224E+09  0.1609E+09  0.2570E+09 
prop  travel  0.1877E+01proj  travel  0.2186E+01 

0.8110E-02  0.2033E+06  0.1060E+04  0.2291E+01  0.2176E+09  0.1572E+09  0.2516E+09 
prop  travel  0.1969E+01proj  travel  0.2291E+01 

0.8210E-02  0.1987E+06  0.1080E+04  0.2398E+O1  0.2130E+09  0.1537E+09  0.2463E+09 
prop  travel  0.2062E+01proj  travel  0.2398E+01 

0.8310E-02  0.1943E+06  0.1100E+04  0.2507E+C1  0.2084E+09  0.1502E+09  0.241 1E+09 
prop  travel  0.2158E+01proj  travel  0.2507E+01 

0.8410E-02  0.1899E+06  0.1119E+04  0.2618E+01  0.2040E+09  0.1469E+09  0.2360E+09 
prop  travel  0.2255E+01proj  travel  0.2618E+01 

0.8510E-02  0.1857E+06  0.1138E+04  0.2731E+01  0.1997E+09  0.1436E+09  0.2310E+09 
prop  travel  0.2354E+01proj  travel  0.2731E+01 

0.8610E-02  0.1817E+06  0.1156E+04  0.2845E+01  0.1954E+09  0.1405E+09  0.2261E+09 
prop  travel  0.2455E+01proj  travel  0.2845E+01 

0.8710E-02  0.1777E+06  0.1174E+04  0.2962E-f01  0.1913E+09  0.1374E+C9  0.2213E+09 
prop  travel  0.2558E+01proj  travel  0.2962E+01 
propellant  1  has  slivered 

0.8810E-02  0.1772E+06  0.1192E+04  0.3080E+01  0.1852E+09  0.1370E+09  0.2116E+09 
prop  travel  0.2662E+01proj  travel  0.3080E+01 

0.8910E-02  0.1773E+06  0.1209E+04  0.3200E+01  0.1814E+09  0.1371E+09  0.2052E+09 
prop  travel  0.2768E+01proj  travel  0.3200E+01 

0.9010E-02  0.1726E+06  0.1227E+04  0.3322E+01  0.1765E+09  0.1335E+09  0.1997E+09 
prop  travel  0.2875E+01proj  travel  0.3322E+01 

0.9110E-02  0.1674E+06  0.1244E+04  0.3446E+01  0.1714E+09  0.1294E+09  0.1938E+09 
prop  travel  0.2982E+01proj  travel  0.3446E+01 

0.9210E-02  0.1621E+06  0.1260E+04  0.3571E+01  0.1662E+09  0.1254E+09  0.1879E+09 
prop  travel  0.3091E+01proj  travel  0.3571E+01 

0.9310E-02  0.1570E+06  0.1276E+04  0.3698E+01  0.1610E+09  0.1214E+09  0.1821E+09 
prop  travel  0.3200E+01proj  travel  0.3698E+01 

0.9410E-02  0.1519E+06  0.1292E+04  0.3826E+01  0.1559E+09  0.1175E+09  0.1763E+09 
prop  travel  0.3310E+01proj  travel  0.3826E+01 
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0.9510E-02  0.1470E+06  0.1307E+04  0.3956E+01  0.1509E+09  0.1137E+09  0.1707E+09 
prop  travel  0.3421E+01proj  travel  0.3956E+01 

0.9610E-02  0.1422E+06  0.1321E+04  0.4087E+01  0.1461E+09  0.1100E+09  0.1652E+09 
prop  travel  0.3533E+01proj  travel  0.4087E+01 

0.9710E-02  0.1376E+06  0.1335E+04  0.4220E+01  0.1414E+09  0.1064E+09  0.1599E+09 
prop  travel  0.3646E+01proj  travel  0.4220E+01 

0.9810E-02  0.1332E+06  0.1349E+04  0.4354E+01  0.1368E+09  0.1030E+09  0.1548E+09 
prop  travel  0.3759E+01proj  travel  0.4354E+01 

0.9910E-02  0.1289E+06  0.1362E+04  0.4490E+01  0.1325E+09  0.9970E+08  0.1498E+09 
prop  travel  0.3873E+01proj  travel  0.4490E+01 

deltat  t  0.998000E-02  intg  t  0.998000E-02 

PMAXMEAN  Pa  0.301051E+09  time  at  PMAXMEAN  sec  0.608000E-02 
PMAXBASE  Pa  0.227560E+09  time  at  PMAXBASE  sec  0.598000E-02 
PMAXBREECH  Pa  0.337155E+09  lime  at  PMAXBREECH  sec  0.618000E-02 
muzzle  VELOCITY  m/s  0.136949E+04  time  of  muzzle  velocity  sec  0.997004E-02 
total  initial  energy  available  J  =  0.483777E+08 

FOR  PROPELLANT  1  MASSFRACT  BURNT  IS  0.966455E+00  at  time  in  sec  0.000000E+00 

total  energy  remaining  in  gas  J=  0.329127E+08 

energy  loss  from  projectile  translation  J=  0.920305E+07 

energy  loss  from  projectile  rotation  J=  0.000000E+00 

energy  lost  to  gas  and  propellant  motion  J=  0.311553E+07 

energy  lost  to  bore  resistance  J=  O.OOOOOOE+OO 

energy  lost  to  recoil  J=  0.1881 16E- 11 

energy  loss  from  heat  transfer  J=  0.152361E+07 

energy  lost  to  air  resistance  J=  O.OOOOOOE+OO 
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APPENDIX  8: 
INPUT  DATA  FOR  XKTC 
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GRADIENT  TESTS  600  IN-  X6315A 
TFFFTTT  110  00 

25  -3  0  3500 

10.0  180.0  0.0001  2.0  0.05  0.01  0.0001  0.0001 

300  100  1100  100  1500  100 

600500120008  0 

0 

529.  14.7  28.896  1.4 

529.0 

JA2  7PF  LOT  RH2-5  0.0  21.3  21.6  .06 

7  .4845  0.020  1.25  7. 


15000. 

1.0 

41754. 

.5 

8000. 

.00030 

1.000 

100000. 

.000300 

1.000 

0.0 

520. 

19826000.  23.0000 

1.2300 

27.00 

0.0 

3.0279 

18.3 

3.0279 

21.3 

2.5 

21.9 

2.5 

23.37 

2.5 

210.56 

2.5 

6.0 

0. 

6.25 

0. 

7.2 

00. 

7.8 

0. 

210.56 

0. 

7.77 

0.0228 

0.7 

21.3 

21.6 

44.0 

0.000 

0.0 

3. 

12. 

30. 

40. 

50. 

80. 

30. 
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No  of 

Copies  Organization 


No  of 

Copies  Organization 


1  Office  of  the  Secretary  of  Defense 
OUSD(A) 

Director,  Live  Fire  Testing 
ATTN:  James  F.  O’Bryon 
Washington,  DC  20301-3110 


1  Director 

US  Army  Aviation  Research 
and  Technology  Activity 
Ames  Research  Center 
Moffett  Field,  CA  94035-1099 


2  Administrator 

Defense  Technical  Info  Center 
ATTN:  DTIC-DDA 
Cameron  Station 
Alexandria,  VA  22304-6145 

1  HQDA  (SARD-TR) 

WASH  DC  20310-0001 

1  Commander 

US  Army  Materiel  Command 
ATTN:  AMCDRA-ST 
5001  Eisenhower  Avenue 

•  4  «  •  «  «  i  AMAAA  AAftl 

Alexandria,  va  zzjjj-wui 


1  Commander 

US  Army  Missile  Command 
ATTN:  AMSMI-RD-CS-R  (DOC) 

Redstone  Arsenal,  AL  35898-5010 

1  Commander 

US  Army  Tank-Automotive  Command 
ATTN:  AMSTA-TSL  (Technical  Library) 
Warren,  MI  48397-5000 

1  Director 

US  Army  TRADOC  Analysis  Command 
ATTN:  ATAA-SL 

While  Sands  Missile  Range,  NM  88002-5502 


1  Commander 

US  Army  Laboratory  Command 
ATTN:  AMSLC-DL 
Adelphi,  MD  20783-1145 


(ciau.  only)  i  Commandant 

US  Army  Infantry  School 
ATTN:  ATSH-CD  (Security  Mgr.) 
Fort  Benning,  GA  31905-5660 


2  Commander 

US  Army,  ARDEC 
ATTN:  SMCAR-IMI-I 
Picatinny  Arsenal,  NJ  07806-5000 


(Undaa.  only)  i  Commandant 

US  Army  Infantry  School 
ATTN:  ATSH-CD-CSO-OR 
Fort  Benning,  GA  31905-5660 


2  Commander 

US  Army,  ARDEC 
ATTN:  SMCAR-TDC 
Picatinny  Arsenal,  NJ  07806-5000 

1  Director 

Benet  Weapons  Laboratory 
US  Army,  ARDEC 
ATTN:  SMCAR-CCB-TL 
Watervliet,  NY  121894050 

1  Commander 

US  Army  Armament,  Munitions 
and  Chemical  Command 
ATTN:  SMCAR-ESP-L 
Rock  Island,  IL  61299-5000 

1  Commander 

US  Army  Aviation  Systems  Command 
ATTN:  AMSAV-DACL 
4300  Goodfellow  Blvd. 

St.  Louis,  MO  63120-1798 


1  Air  Force  Armament  Laboratory 
ATTN:  AFATL/DLODL 
Eglin  AFB,  FL  32542-5000 

Aberdeen  Proving  Ground 

2  Dir,  USAMSAA 

ATTN:  AMXSY-D 

AMXSY-MP,  H.  Cohen 

1  Cdr,  USATECOM 

ATTN:  AMSTE-TD 

3  Cdr,  CRDEC,  AMCCOM 

ATTN:  SMCCR-RSP-A 
SMCCR-MU 
SMCCR-MSI 

1  Dir,  VLAMO 

ATTN:  AMSLC-VL-D 
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No.  of 

No.  of 

Codes 

Organization 

Copies 

Organization 

1 

Commander 

USA  Concepts  Analysis  Agency 

ATTN:  D.  Hardison 

8120  Woodmont  Avenue 

Bethesda,  MD  20014-2797 

3 

PEO-Armaments 

Project  Manger 

Tank  Main  Armament  Systems 

ATTN:  AMCPM-TMA,  K.  Russell 
AMCPM-TM  A- 105 
AMCPM-TMA-120 

1 

C.IA 

OIR/DB/Standard 

Picatinny  Arsenal,  NJ  07806-5000 

Washington,  DC  20505 

1 

Commander 

Armament  RD&E  Center 

l 

US  Army  Ballistic  Missile 

Defense  Systems  Command 

Advanced  Technology  Center 

P.O.  Box  1500 

US  Army  AMCCOM 

ATTN:  SMCAR-AEE 

Picatinny  Arsenal,  NJ  07806-5000 

Huntsville,  AL  35807-3801 

8 

Commander 

Armament  RD&E  Center 

l 

Chairman 

DoD  Explosives  Safety  Board 

Room  856-C 

Hoffman  Bldg.  1 

2461  Eisenhower  Avenue 

Alexandria,  VA  22331-0600 

US  Army  AMCCOM 

ATTN:  SMCAR-AEE-B 

A.  Beardell 

B.  Biodman 

D.  Downs 

S.  Einstein 

S.  Westley 

l 

Commander 

US  Army  Material  Command 

ATTN:  AMCPM-GCM-WF 

5001  Eisenhower  Avenue 

Alexandria,  VA  22333-5001 

2 

S.  Bernstein 

C.  Roller 

J.  Rutkowski 

Picatinny  Arsenal,  NJ  07806-5000 

Commander 

l 

Commander 

US  Army  Material  Command 

ATTN;  AMCDE-DW 

5001  Eisenhower  Avenue 

Alexandria,  VA  22333-5001 

2 

US  Army  ARDEC 

ATTN:  SMCAR-AES,  S.  Kaplowitz; 

D.  Spring 

Picatinny  Arsenal,  NJ  07806-5000 

Commander 

4 

Project  Manager 

Autonomous  Precision-Guided 

Munition  (APGM) 

Armament  RD&E  Center 

US  Army  AMCCOM 

ATTN:  AMCPM-CW 

Armament  RD&E  Center 

US  Army  AMCCOM 

ATTN:  SMCAR-HFM,  E.  Barrieres 
SMCAR-CCH-V,  C.  Mandala 
Picatinny  Arsenal,  NJ  07806-5000 

2 

AMCPM-CWW 

AMCPM-CWS,  M.  Fisette 
AMCPM-CWA-S,  R.  DeKleine 

Pica  tinny  Arsenal,  NJ  07806-5000 

Project  Manager 

1 

Commander 

Armament  RD&E  Center 

US  Army  AMCCOM 

ATTN:  SMCAR-FSA-T,  M.  Salsbury 
Picatinny  Arsenal,  NJ  07806-5000 

Production  Base  Modernization  Agency 
ATTN:  AMSMC-PBM,  A.  Siklosi 

AMSMC-PBM-E,  L.  Laibson 
Picatinny  Arsenal,  NJ  07806-5000 

1 

Commander,  USACECOM 

CECOM  R&D  Technical  Library 
ATTN:  ASQNC-ELC-I-T,  Myer  Center 
Fort  Monmouth,  NJ  07703-5301 
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No.  of 

Copies  Organization 

1  Commander 

US  Army  Harry  Diamond  Laboratories 
ATTN:  SLCHD-TA-L 
2800  Powder  Mill  Rd 
Adelphi,  MD  20783-1145 

1  Commandant 

US  Army  Aviation  School 
ATTN:  Aviation  Agency 
Fort  Rucker,  AL  36360 

1  Project  Manager 

US  Army  Tank-Automotive  Command 
Improved  TOW  Vehicle 
ATTN:  AMCPM-ITV 
Warren,  MI  48397-5000 

2  Program  Manager 

US  Army  Tank-Automotive  Command 
ATTN:  AMCPM-ABMS,  T.  Dean 
Warren,  MI  48092-2498 

1  Project  Manager 

US  Army  Tank-Automotive  Command 
Fighting  Vehicle  Systems 
ATTN:  AMCPM-BFVS 
Warren,  MI  48092-2498 

1  President 

US  Army  Armor  and  Engineer  Board 
ATTN:  ATZK-AD-S 
Fort  Knox,  KY  40121-5200 

1  Project  Manager 

US  Army  Tank-Automotive  Command 
M-60  Tank  Development 
ATTN:  AMCPM-ABMS  (PM  Abrams) 
Wanen,  MI  48092-2498 

1  Commander 

US  Army  Training  and  Doctrine  Command 
ATTN:  ATCD-MA,  MAJ  William® 

Fort  Monroe.  VA  23651 

2  Director 

US  Army  Materials  Technology  Laboratory 
ATTN:  SLCMT-ATL 
Watertown,  MA  02172-0001 

1  Commander 

US  Army  Research  Office 
ATTN:  Technical  Library 
P.  O.  Box  12211 

Research  Triangle  Park,  NC  27709-2211 


No.  of 

Copies  Organization 

1  Commander 

US  Army  Belvoir  Research  and 
Development  Center 
ATTN:  STRBE-WC 
Fort  Belvoir,  VA  22060-5006 

1  Director 

US  Army  TRAC-Ft  Lee 
ATTN:  ATRC-L,  (Mr.  Cameron) 
Fort  L oe,  VA  23801-6140 

1  Commandant 

US  Army  Command  and  General 
Staff  College 

Fort  Leavenworth,  KS  66027 

1  Commandant 

US  Army  Special  Warfare  School 
ATTN:  Rev  and  Tng  Lit  Div 
Fort  Bragg,  NC  28307 

3  Commander 

Radford  Army  Ammunition  Plant 
ATTN:  SMCAR-QA/HI  LIB 
Radford,  VA  24141-0298 

1  Commander 

US  Army  Foreign  Science  and 
Technology  Center 
ATIN:  AMXST-MC-3 
220  Seventh  Street,  NE 
Charlottesville,  VA  22901-5396 

2  Commander 

Naval  Sea  Systems  Command 
ATTN:  SEA62P 
SEA  64 

Washington,  DC  20362-5101 

1  Commander 

Naval  Air  Systems  Command 
ATTN:  AIR-954-Technical  Library 
Washington,  DC  20360 

1  Assistant  Secretary  of  the  Navy 

(R,  E,  and  S) 

ATIN:  R.  Reichenbach 
Room  5E787 
Pentagon  Bldg 
Washington,  DC  20375 
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No.  of 

Copies  Organization 

1  Naval  Research  Laboratory 

Technical  Library 
Washington.  DC  20375 

1  Commandant 

US  Army  Command  and  General 
Staff  College 

Fort  Leavenworth,  KS  66027 

2  Commandant 

US  Army  Field  Artillery  Cotter 
and  School 

ATTN:  ATSF-CO-MW,  B.  Willis 
Ft.  Sill.  OK  73503-5600 

1  Office  of  Naval  Research 
ATTN:  Code  473,  R.  S.  Miller 
800  N.  Quincy  Street 
Arlington,  VA  22217-9999 

3  Commandant 

US  Army  Armor  School 

ATTN:  ATZK-CD-MS,  M.  Falkoviteh 

Armor  Agency 

Fort  Knox,  KY  40121-5215 

2  Commander 

US  Naval  Surface  Warfare  Center 
ATTN:  J.  P.  Consaga 
C.  Gotzmer 

Indian  Head,  MD  20640-5000 

4  Commander 

Naval  Surface  Warfare  Cotter 
ATTN:  Code  240,  S.  Jacobs 
Code  730 

Code  R-13,  K.  Kim 
R.  Bemecker 

Silver  Spring,  MD  20903-5000 

2  Commanding  Officer 

Naval  Underwater  Systems  Center 
ATTN:  Code  5B331,  R.  S.  Lazar 
Technical  Library 
Newport,  RI  02840 

5  Commander 

Naval  Surface  Warfare  Center 
ATTN:  Code  G33,  J.  L.  East 
W.  Burrell 
J.  Johndrow 
Code  G23,  D.  McClure 
Code  DX-21  Technical  Library 
Dahlgren,  VA  22448-5000 


No.  of 

Copies  Organization 

3  Commander 

Naval  Weapons  Center 
ATTN:  Code  388,  C.  F.  Price 
Code  3895,  T.  Parr 
Information  Science  Division 
China  Lake,  CA  93  *'55-6001 

1  Program  Manager 

AFOSR 

Directorate  of  Aerospace  Sciences 
ATTN:  L.  H.  Caveny 
Bolling  AFB 

Washington,  DC  20332-0001 

5  Commander 

Naval  Ordnance  Station 
ATTN:  L.  Torreyson 
T.  C.  Smith 
D.  Brooks 
W.  Vienna 
Technical  Library 
Indian  Head,  MD  20640-5000 

1  AL/TSTL  (Technical  Library) 

ATTN:  J.  Lamb 

Edwards  AFB,  CA  93523-5000 

1  AFSC/SDOA 

Andrews  AFB,  MD  20334 

1  AFATL/DLY  V 

Eglin  AFB,  FL  32542-5000 

1  AFATL/DLXP 

Eglin  AFB,  FL  32542-5000 

1  AFATL/DUE 

Eglin  AFB,  FL  32542-5000 

1  NASA/Lyndon  B.  Johnson  Space  Center 

ATTN:  NHS-22  Library  Section 
Houston,  TX  77054 

1  AFELM,  The  Rand  Corporation 

ATTN:  Library  D 
1700  Main  Street 
Santa  Monica,  CA  90401-3297 

3  AAI  Corporation 

ATTN:  J.  Herbert 
J.  Frankie 
D.  Cleveland 
P.  O.  Box  126 

Hunt  Valley,  MD  21030-0126 
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No.  of 

Copies  Organization 

1  Aerojet  Ordnance  Company 

ATTN:  D.  Thatcher 
2521  Michelle  Drive 
Tustin,  CA  92680-7014 

• 

1  Aerojet  Solid  Propulsion  Company 

ATTN:  P.  Micheli 
Sarcramento,  CA  96813 

1  Atlantic  Research  Corporation 

ATTN:  M.  King 
5390  Cherokee  Avenue 
Alexandria,  VA  22312-2302 

4  AL/LSCF 

ATTN:  J.  Levine 
L.  Quinn 
D.  Williams 
T.  Edwards 

Edwards  AFB,  CA  93523-5000 

1  AVCO  Everett  Research  Laboratory 
ATTN:  D.  Stickler 

2385  Revere  Beach  Parkway 
Everett,  MA  02149-5936 

2  Cal  span  Corporation 
ATTN:  C.  Murphy 
P.  O.  Box  400 
Buffalo,  NY  14225-0400 

1  General  Electric  Company 

Armament  Systems  Department 
ATTN:  M.  J.  Bulman 
128  Lakeside  Avenue 
Burlington,  VT  054014985 

1  IITRI 

ATTN:  M.  J.  Klein 
10  W.  35th  Street 
Chicago,  IL  60616-3799 

1  Hercules,  Inc. 

Allegheny  Ballistics  Laboratory 
ATTN:  William  B.  Walkup 
P.  O.  Box  210 
Rocket  Center,  WV  26726 

1  Hercules,  Inc. 

Radford  Army  Ammunition  Plant 
ATTN:  J.  Pierce 
Radford,  VA  24141-0299 


No.  of 

Copies  Organization 

3  Lawrence  Livermore  National 

Laboratory 

ATTN:  L-355,  A.  Buckingham 
M.  Finger 

L-324,  M.  Constantino 
P.  O.  Box  808 
Livermore  CA  94550-0622 

1  Olin  Corporation 

Badger  Army  Ammunition  Plant 
ATTN:  R.  J.  Thiede 
Baraboo,  WI  53913 

1  Olin  Corporation 

Smokeless  Powder  Operation 
ATTN:  D.  C.  Mann 
P.  O.  Box  222 
Si.  Marks,  FL  32355-0222 

1  Paul  Gough  Associates,  Inc. 

ATTN:  P.  S.  Gough 
P.  O.  Box  164 
1048  South  Street 
Portsmouth,  NH  03801-1614 

1  Physics  International  Company 

ATTN:  Library,  H.  Wayne  Wampler 

2700  Merced  Street 

San  Leandro,  CA  98457-5602 

1  Princeton  Combustion  Research 

Laboratory,  Inc. 

ATTN:  M.  Summerfield 

475  US  Highway  One 

Monmouth  Junction,  NJ  08852-9650 

2  Rockwell  International 
Rocketdyne  Division 
ATTN:  BA08,  J.  £.  Flanagan 

J.  Gray 

6633  Canoga  Avenue 
Canoga  Park,  CA  91303-2703 

3  Thiokol  Corporation 
Huntsville  Division 
ATTN:  D.  Flanigan 

Dr.  John  Deur 
Technical  Library 
Huntsville,  AL  35807 
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No.  of 

No.  of 

Copies 

Organization 

Copies 

Organization 

2 

Thiokol  Corporation 

Elkton  Division 

ATTN:  R.  Biddle 

Technical  Library 

P.  0.  Box  241 

1 

University  of  Minnesota 

Department  of  Mechanical  Engineering 
ATTN:  E.  Fletcher 

Minneapolis,  MN  55414-3368 

Elkton,  MD  21921-0241 

1 

Case  Western  Reserve  University 
Division  of  Aerospace  Sciences 

1 

Veritay  Technology,  Inc. 

ATTN:  E.  Fisher 

4845  Millersport  Highway 

ATTN:  J.  Tien 

Cleveland,  OH  44135 

East  Amherst,  NY  14501-0305 

3 

Georgia  Institute  of  Technology 

School  of  Aerospace  Engineering 

1 

Universal  Propulsion  Company 

ATTN:  H.  J.  McSpadden 

Black  Canyon  Stage  1 

Box  1140 

Phoenix,  AZ  84029 

1 

ATTN:  B.T.  Zim 

E.  Price 

W.C.  Strahle 

Atlanta,  GA  30332 

Institute  of  Gas  Technology 

1 

Battelle  Memorial  Institute 

ATTN:  Technical  Library 

505  King  Avenue 

Columbus,  OH  43201-2693 

1 

ATTN:  D.  Gidaspow 

3424  S.  State  Street 

Chicago.  IL  60616-3896 

Johns  Hopkins  University 

1 

Brigham  Young  University 

Department  rtf  Chemical  Engineering 

ATTN:  M.  Beckstead 

Provo,  UT  84601 

Applied  Physics  Laboratory 

Chemical  Propulsion 

Information  Agency 

ATTN:  T.  Christian 

Johns  Hopkins  Road 

1 

California  Institute  of  Technology 

204  Karman  Laboratory 

Laurel,  MD  20707-0690 

1 

Main  Stop  301-46 

ATTN:  F£.C.  Culick 

1201  E.  California  Street 

Pasadena,  CA  91109 

California  Institute  of  Technology 

1 

Massachusetts  Institute  of  Technology 
Department  of  Mechanical  Engineering 
ATTN:  T.  Toong 

77  Massachusetts  Avenue 

Cambridge,  MA  021394307 

Jet  Propulsion  Laboratory 

ATTN:  L.  D.  Strand,  MS  512/102 

4800  Oak  Grove  Drive 

Pasadena,  CA  91109-8099 

1 

Pennsylvania  State  University 

Applied  Research  Laboratory 

ATTN:  G.M.  Faeth 

University  Park,  PA  16802-7501 

1 

University  of  Illinois 

Department  rtf  Mechanical/Industrial 
Engineering 

ATTN:  H.  Krier 

144  MEB;  1206  N.  Great  Street 

1 

Pennsylvania  State  University 
Department  of  Mechanical  Engineering 
ATTN:  K.  Kuo 

University  Park,  PA  16802-7501 

Urbana,  IL  61801-2978 

1 

Purdue  University 

School  of  Mechanical  Engineering 

1 

University  of  Massachusetts 

Department  of  Mechanical  Engineering 
ATTN:  K.  Jakus 

Amherst,  MA  01002-0014 

ATTN:  J.  R.  Osborn 

TSPC  Chaffee  Hall 

West  Lafayette,  IN  47907-1199 
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Copies  Organization 

1  SRI  International 

Propulsion  Sciences  Division 
ATTN:  Technical  Library 
333  Ravenwood  Avenue 
Menlo  Park,  CA  94025-3493 

1  Rensselaer  Ploytechnic  Institute 

Department  of  Mathematics 
Troy,  NY  12181 

2  Director 

Los  Alamos  Scientific  Laboratory 
ATTN:  T3,D.  Butler 

M.  Division,  B.  Craig 
P.  O.  Box  1663 
Los  Alamos,  NM  87544 

1  General  Applied  Sciences  Laboratory 

ATTN:  J.  Erdos 

77  Raynor  Avenue 

Ronkonkama,  NY  11779-6649 

1  Battelle  PNL 

ATTN:  Mr.  Mark  Gamich 
P.  O.  Box  999 
Richland,  WA  99352 

1  Stevens  Institute  of  Technology 

Davidson  Laboratory 
ATTN:  R.  McAlevy,  III 
Castle  Point  Station 
Hoboken,  NJ  07030-5907 

1  Rutgers  University 

Department  of  Mechanical  and 
Aerospace  Engineering 
ATTN:  S.  Temkin 
University  Heights  Campus 
New  Brunswick,  NJ  08903 

1  University  of  Southern  California 

Mechanical  Engineering  Department 
ATTN:  0HE200,  M.  Gerstein 

Los  Angeles,  CA  90089-5199 

2  University  of  Utah 

Department  of  Chemical  Engineering 
ATTN:  A.  Baer 

G.  Flandro 

Salt  Lake  City,  UT  84112-1194 


No.  of 

Copies  Organization 

1  Washington  State  University 

Department  of  Mechanical  Engineering 
ATTN:  C.  T.  Crowe 
Pullman,  WA  99163-5201 

1  Honeywell,  Inc. 

ATTN:  R.  E.  Tompkins 
MN38-3300 

10400  Yellow  Circle  Drive 
Minnetonka,  MN  55343 

1  Science  Applications,  Inc. 

ATTN:  R.  B.  Edelman 
23146  Cumorah  Crest  Drive 
Woodland  Hills,  CA  91364-3710 

Aberdeen  Proving  Ground 

Cdr,  CSTA 

ATTN:  STECS-LI,  R.  Hendriclcsen 
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